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Executive  Summary 

This  document  is  the  summary  of  the  final  outcome  of  the  research  contract  w7714- 
1 15195/A  entitled  Research  &  development  for  Long  Term  Evolution  (LTE)  wireless 
location  based  on  synthetic  array.  The  report  is  a  compilation  of  the  research  activities 
and  outcomes  for  the  period  of  January  2012  to  March  31  2014. 

The  final  report  consists  of  two  parts  whereas  this  is  part  B.  This  document  is  primarily 
based  on  the  PhD  dissertation  of  Dr.  Neda  Moazen  who  did  her  doctoral  research  as 
part  of  this  contract. 
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Acronyms 


AN 

Access  node  (eNodeB) 

AP 

Anchor  point 

ANan 

Access  node  that  is  an  anchor  node 

ANfp 

Access  node  that  is  feature  point 

AWGN 

Additive  White  Gaussian  Noise 

AOA 

Angle  of  Arrival 

BCRLB 

Bayesian  Cramer-Rao  Lower  Bound 

BF 

Bayesian  Filter 

BFI 

Bayesian  Fisher  Information 

BFIM 

Bayesian  Fisher  Information  Matrix 

CDMA 

Code  Division  Multiple  Access 

CML 

Concurrent  Mapping  and  Localization 

CRLB 

Cramer-Rao  Lower  Bound 

CV 

Computer  Vision 

EKF 

Extended  Kalman  Filter 

FastSLAM 

Factorized  Solution  for  Simultaneous  Localization  and  Mapping 

FI 

Fisher  Information 

FIM 

Fisher  Information  Matrix 

FP 

Feature  Point 

GNSS 

Global  Navigation  Satellite  Systems 

GPS 

Global  Positioning  System 

i.i.d 

Independent  and  Identically  Distributed 

IMM 

Interacting  Multiple  Model 

IMU 

Inertial  Measurement  Unit 
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IEEE 

Institute  of  Electrical  and  Electronics  Engineers 

KF 

Kalman  Filter 

KL 

Kullback-Leibler 

LAMBDA 

Least-squares  Ambiguity  Decorrelation  Adjustment 

LJG 

Linear  and  Jointly  Gaussian 

LML 

Local  Maximum  Likelihood 

LOS 

Line  of  Sight 

LTE 

Long  Term  Evolution 

MAP 

Maximum  A  Posterior 

MC 

Monte  Carlo 

M-CRLB 

Modified  Cramer-Rao  Lower  Bound 

MEMS 

Micro  Electro-Mechanical  Systems 

Ml 

Mutual  Information 

MLE 

Maximum  Likelihood  Estimator 

MM 

Multiple  Model 

MMSE 

Minimum  Mean  Square  Error 

MN 

Mobile  Node 

NEES 

Normalized  Estimation  Error  Squared 

OCXO 

Oven-Controlled  Crystal  Oscillator 

OWLS 

Opportunistic  Wireless  Localization  System 

PDF 

Probability  Density  Function 

PDOA 

Phase  Difference  of  Arrival 

PEB 

Position  Error  Bound 

PF 

Particle  Filter 

PRN 

Pseudo-Random  Noise 
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POA 

Phase  of  Arrival 

RB 

Rao-Blackwellized 

RBPF 

Rao-Blackwellized  Particle  Filter 

RFID 

Radio  Frequency  Identification 

RIM 

Research  in  Motion 

RTOF 

Round-trip  Time  of  Flight 

SLAM 

Simultaneous  Localization  and  Mapping 

SNR 

Signal-to-Noise  Ratio 

SS 

Signal  Strength 

TCXO 

Temperature  Compensated  Crystal  Oscillator 

TDOA 

Time  Difference  of  Arrival 

TOA 

Time  of  Arrival 

UWB 

Ultra  Wideband 

WiFi 

Wireless  Fidelity 

WLAN 

Wireless  Local  Access  Network 

Variables  used 


variable 

description 

bel(x0:t) 

Posterior  pdf  of  the  particles  over  the  complete  time  interval  from  0  to  tTs 

b, 

range  offset 

Bs 

signal  bandwidth 

c 

propagation  velocity 

C1  xm 

1  xm  row  vector  whose  elements  are  zero,  except  the  i -th  element  which  is 

equal  to  one. 
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/ 

carrier  frequency 

fb 

clock  frequency  bias 

1 

clock  frequency  drift 

m 

clock  random  frequency  error 

I(X,Y) 

mutual  information  between  Y  and  X 

J 

Fisher  information  matrix  in  general 

^ tot 

total  Bayesian  Fisher  information  matrix 

J2 

measurement  information  matrix 

JP 

apriori  information  matrix 

Kan 

vector  of  direction 

m 

stacked  vector  of  AN  locations 

m 

state  vector  describing  the  location  of  the  /th  ANs 

stacked  vector  of  APs  location 

mfP 

stacked  vector  of  FPs  location 

{mt*’mUy} 

2D  location  variables  of  an  AN 

number  of  ANs 

Nap 

number  of  APs 

Nfp 

number  of  FPs 

Nd 

dimension  of  dynamic  variable  in  state  vector 

Ns 

dimension  of  stationary  variable  in  state  vector 

NP 

number  of  particles 

V, 

MN  location  vector 

Pi, 

history  of  MN  locations 

q, 

state  vector  at  time  step  t 
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dynamic  portion  of  state  vector  at  time  step  t 

q. 

stationary  portion  of  state  vector 

a 

covariance  matrix  of  dynamic  variables  update  process 

0, 

motion  process  covariance  matrix 

r(t) 

geometric  range  between  the  AN  and  MN 

s(t) 

bandpass  signal 

~s(t) 

lowpass  signal 

U, 

control  vector 

u,:, 

history  of  update  control  inputs 

V, 

motion  update  process  noise 

wi,< 

measurement  noise 

Zr 

measurement  vector  received  at  time  step  t 

Zi,t 

observation  vector  from  /- th  AN  at  the  time  t 

history  of  observation  from  /-th  AN 

SiJ 

uk+ 1 

{k  +  \)y.(k  + 1)  dimensional  matrix  whose  elements  are  all  zero  except  at  the  / - 

th  row  and  the  j  -th  column  which  is  one 

fit MN 

MN  clock  bias 

fit  AN 

AN  clock  bias 

fi^MN 

MN  clock  rate  drift 

fit  AN 

AN  clock  rate  drift 

MO 

carrier  phase  variation 

A'8 

a 

VaV/ ,  second  order  derivative  operator 

carrier  phase  measurement  noise 

n 

normalizing  constant 
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A 

carrier  wavelength 

N(b,B) 

Gaussian  process  with  the  mean  vector  b  and  covariance  matrix  B 

n 

overall  transition  matrix  of  sight  states 

n, 

transition  matrix  of  the  /'- th  AN’s  sight  state 

°l(Qb) 

range  offset  variance  (covariance  matrix) 

% 

carrier  phase  of  AN  transmitter  at  time  of  transmission 

<p(t) 

partial  carrier  phase  cycle  measurement 

x2 

chi-square  distribution 

dim(.) 

dimension  function 

Re{.} 

real  part 

V. 

Jacobian  function 
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Chapter  1  Introduction 


In  recent  years,  new  developments  in  wireless  systems  have  combined  mobility 
features  with  personal  communication  devices,  such  as  laptops  and  cell  phones.  However, 
advances  in  embedded  sensing  and  high  speed  processing  have  enabled  location-aware 
services  built  into  handset  devices.  The  most  popular  developed  positioning  systems  used  by 
handset  devices  are  based  on  the  Global  Positioning  System  (GPS)  and  cellular  networks  which 
are  widely  accessible  for  navigation  in  outdoor  environments. 

Each  GPS  satellite  transmits  a  narrowband  signal  modulated  with  a  unique  pseudo¬ 
random  noise  (PRN)  sequence  code  that  is  known  to  the  receiver.  If  the  transmitter  and 
receivers  are  synchronized,  the  corresponding  receiver  can  correlate  the  received  signal  with  a 
known  PRN  sequence  to  estimate  the  propagation  delay  between  the  transmitter  and  receiver 
antennas.  The  typical  GPS  navigation  solution  requires  a  minimum  of  4  satellites  that  must  be 
simultaneously  visible  to  estimate  the  3D  location  of  the  receiver  and  the  time  offset  between 
receiver  and  transmitter  clocks.  However,  in  indoor  and  dense  urban  environments,  satellite 
visibility  is  problematic. 

On  the  other  hand,  the  position  accuracy  obtained  by  existing  cellular-based  methods 
using  cell-1  D  or  enhanced  observed  time  difference  of  arrival  (E-OTDOA)  is  generally  low  and 
insufficient  for  indoor  environments.  To  this  end,  indoor  localization  systems  rely  on  other 
technologies  such  as  UWB,  RFID,  Bluetooth  (IEEE  802.15),  ultrasonic  badges,  and  computer 
vision  technology.  However,  these  approaches  require  developing  a  signaling  system  and/or 
installing  the  network  infrastructure  which  make  them  time-consuming  and  costly  approaches 
for  indoor  localization. 

Particularly,  the  project  of  “Google  Maps”  extension  inside  buildings  has  sparked  a 
technology  race  among  all  smartphone  manufacturers  such  as  Apple,  Nokia  ,  Samsung,  and 
RIM;  and  ASIC  manufactures  such  as  Qualcomm  and  Broadcom.  This  is  augmented  by  the 
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concept  of  using  any  opportunistic  wireless  signal  that  may  be  present.  In  this  work,  the 
opportunistic  wireless  localization  system  (OWLS)  is  developed  where  the  wireless  signals 
utilized  for  positioning  are  received  from  opportunistic  and  poorly  defined  sources  that  perhaps 
are  uncooperative  with  the  mobile  handset  device.  This  means  that  the  mobile  node  (MN)  is  not 
registered  with  the  network  and  merely  exploits  received  signals  as  signals  of  opportunity. 

This  research  has  focused  on  the  OWLS  using  existing  wireless  network  infrastructure 
such  as  locally  generated  WiFi,  WLAN  and  4G  LTE  wireless  signals  as  an  inexpensive  solution 
for  indoor  and  dense  urban  environments.  In  addition  to  cost-effectiveness,  opportunistic 
wireless  localization  offers  scalability,  not  only  in  the  cost  of  required  network  infrastructure,  but 
also  in  the  number  of  mobile  devices/nodes  (MNs)  subscribing  to  positioning  services,  since 
each  mobile  device  is  responsible  for  its  sensing  and  processing. 

If  the  signal  strength  map  is  known,  WiFi  fingerprinting  methods  using  signal  of 
opportunity  can  achieve  high  localization  accuracy  in  indoor  environments  [i],  [ii],  [i]  has 
obtained  accuracy  of  0.25  m  using  Monte  Carlo  localization  based  on  a  spatial  discretized  map 
of  signal  strength,  combined  with  contact  sensing.  In  [ii],  the  spatial  discretized  map  of  signal 
strength  is  fused  with  a  low-cost  image  sensor  to  obtain  a  localization  accuracy  of  3  m. 
However,  fingerprinting  methods  rely  on  expensive  and  time-consuming  training  phases  to 
obtain  the  signal  strength  map.  These  methods  require  that  the  area  must  already  have  been 
surveyed  and  this  makes  them  inefficient  approaches  for  real  time  OWLS. 

An  ideal  OWLS  does  not  rely  on  prior  knowledge  of  a  signal  map,  but  builds  its  own 
database  during  operation.  Nevertheless,  in  absence  of  pre-existing  information  about  the 
environment  and  the  signal  strength  map,  an  OWLS  must  deal  with  reception  of  a  multitude  of 
signals  with  a  large  number  of  unknown  channel  and  source  parameters.  Opportunistic  wireless 
localization  introduces  the  simultaneous  localization  and  mapping  problem  (SLAM)  as  a 
systematic  approach  to  not  only  localize  the  MN  but  also  to  estimate  and  track  the  unknown 
parameters  due  to  signal  reception  from  unregistered  networks,  including  the  AN  locations  [iii], 
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[ivj.  SLAM,  also  known  as  Concurrent  Mapping  and  localization  (CML),  is  one  of  the  most 
fundamental  processing  algorithms  in  robotics  that  has  especially  attracted  immense  attention  in 
automation  and  control  of  the  Mars  Exploration  Rover  [v],  [vi]. 

An  OWLS  must  be  solved  for  a  large  number  of  constraints  applied  while  the  MN  moves 
in  an  unknown  trajectory  and  receives  signal  from  unknown  sources.  The  most  significant 
concept  from  SLAM  that  can  be  applied  to  OWLS  is  the  introduction  of  a  systematic  way  of 
incorporating  all  of  the  disparate  nonlinear  constraints  of  the  MN  dynamic  motion  and 
observables.  The  constraints  define  correlation  between  environment  map  parameters  and  MN 
locations  estimates.  Correlation  always  increases  in  time  as  more  constraints  are  defined  while 
the  MN  moves  and/or  receives  new  observations.  The  SLAM  behavior  in  an  OWLS  is 
analogous  to  a  network  of  springs  where  each  constraint  corresponds  to  a  spring  between 
unknowns  [iii],  as  illustrated  in  Figure  0-1.  An  observation  from  an  access  node,  at  each  time 
step,  is  analogous  with  a  displacement  in  a  spring  network,  and  its  correlation  effect  on  its 
neighbors  is  proportional  to  their  distances  to  other  access  nodes,  i.e.,  the  nearer  they  are,  the 
greater  is  the  effect.  The  minimum  energy  solution  for  this  nonlinear  network  of  springs  self- 
organize  the  map  parameters  and  the  MN  trajectory.  As  the  MN  moves  through  the 
environments  and  receives  more  observation,  the  spring  becomes  stiffer  and  fixes  the  map 
parameters. 
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p 


Figure  0-1  SLAM  as  a  network  of  springs  problem 


Additionally,  ANs’  clocks  may  not  necessarily  be  synchronized  such  as  FDD  LTE 
networks.  The  biased  ranging  issue  is  aggravated  in  a  non-stationary  multipath  propagation 
environment  where  the  signal  is  subject  to  spatial  fading  and  signal  attenuation  due  to 
shadowing  and  line-of  sight  (LOS)  blockage  [vii],  [viii]. 

Furthermore,  unlike  traditional  wireless  localization  systems,  the  wireless  source 
locations  are  not  necessarily  unknown  when  the  MN  tends  to  estimate  its  location  based  on  the 
signal  receptions  from  an  unregistered  wireless  network.  Hence,  the  AN  location  parameters 
need  to  be  estimated,  along  with  the  MN  position  states,  which  degrade  the  MN  position 
estimates  due  to  the  increase  in  the  number  of  unknown  variables,  while  observables  are 
unchanged.  Tracking  algorithms  must  be  transformed  to  a  solution  for  a  SLAM  problem  which 
the  MN  localization  is  dependent  on  the  quality  of  estimation  of  AN  locations  [iii],  [iv]. 

The  SLAM-based  OWLS  must  deal  with  increase  of  unknown  random  variables  such  as 
the  MN  trajectory  and  the  range  error  as  the  MN  moves  and  receives  new  observations. 
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Consequently,  the  large  set  of  unknown  variables  will  be  unobservable  since  the  number  of 
unknowns  increases  faster  than  the  number  of  known  equations.  It  will  be  shown  the  OWLS 
algorithm  converges  and  is  practically  robust  if  two  rather  benign  assumptions  can  be  made: 

1-  The  AN’s  are  stationary  in  terms  of  physical  position. 

2-  The  MN  undergoes  a  smooth  albeit  random  trajectory. 

The  assumption  of  smooth  trajectory  is  significant  because  there  is  high  correlation 
between  the  positions  that  can  be  used  to  filter  out  noisy  and  false  location  estimates  due  to  the 
multipath  and  clock  instability  [ixj.  The  Bayesian-based  filter  solution  is  proposed  to  improve  the 
MN  position  estimation  as  it  can  take  advantage  of  correlated  information  from  previous 
observations.  However,  localization  performance  relies  on  the  knowledge  of  the  network  map, 
the  MN  trajectory  and  quality  of  measurements. 

Second  stage  involves  obtaining  an  error  metric  to  assess  the  system  performance,  as 
addressed  in  Chapter  3.  In  estimation  theory,  minimum  mean  square  error  (MMSE),  calculated 
through  the  covariance  matrix  of  posterior  PDF,  is  generally  considered  as  a  measure  of 
estimation  accuracy.  However,  information  theory  can  provide  a  deep  analytical  perspective  to 
assess  the  information  of  observables  regarding  state  space.  One  of  the  main  components  of 
this  research  is  dedicated  to  determining  the  envelope  of  reasonably  robust  operability  of  the 
SLAM  algorithms  for  opportunistic  wireless  localization.  The  Bayesian  Fisher  information  matrix 
is  a  key  to  the  SLAM  problem,  representing  a  quantitative  summary  of  all  available  information 
in  one  matrix.  For  the  SLAM-based  OWLS,  the  Bayesian  Fisher  information  matrix  (BFIM)  is 
derived  as  a  quantitative  measure  of  all  available  information  from  observables,  the  MN  motion 
process  and  prior  knowledge  regarding  the  unknowns.  BFIM  and  its  variants  are  utilized  to 
analyze  the  proposed  SLAM-based  solutions  for  OWLS  performance  in  terms  of  efficiency, 
observability  and  consistency.  Even  though,  the  estimation  relation  to  quantifiable  Fisher 
information  is  only  valid  for  a  jointly  linear  Gaussian  problem;  as  long  as  the  nonlinearity  is  mild 
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and  the  measurement  and  state  update  noise  are  approximately  Gaussian,  it  can  still  be 
regarded  as  a  fundamental  limit  of  accuracy  [x]. 

To  understand  how  much  information  each  unknown  variable  contributes  to  the  problem, 
SLAM-based  OWLS  is  studied  by  four  system  models  according  to  ranging  measurements, 
varying  from  a  simple  synchronized  scenario  with  stationary  known  propagation  to  the  worse 
case  with  unsynchronized  reception  from  unregistered  and  unsynchronized  network  in  a  non¬ 
stationary  mixed  LOS/NLOS  propagation  scenario. 

Input-output  mutual  information  (Ml)  is  another  interesting  measure  of  information  which 
quantifies  the  mutual  dependence  between  input  and  output  variables.  Guo  et  al. ,  2005, 
revealed  an  insightful  relationship  between  mutual  information  and  MMSE  in  the  Gaussian 
channel,  regardless  of  the  input  statistics,  given  by  [xi]: 

dliX-Y)  1  .  .  0-1 

v  ’  =  -MMSE  ( p) 
dp  2 

where  l(X;Y )  is  the  mutual  information  between  input  X  ,  and  output  Y  ,  in  units  of  nats;  and 
p  is  the  signal-to-noise  ratio.  Guo’s  theorem  illustrates  that  mutual  information  as  well  as 
Fisher  information  are  powerful  metrics  which  can  lead  to  different  insight  into  the  localization 
performance.  In  this  regard,  the  application  of  Ml  for  localization  problems  with  special  focus  on 
our  SLAM-based  problem  is  investigated  in  Chapter  3. 

In  a  broad  perspective,  the  SLAM-based  OWLS  can  be  viewed  as  a  dynamic  problem 
varying  from  mapping-only  to  tracking-only,  or  simultaneous  localization  and  mapping, 
depending  on  the  available  knowledge  from  source  locations  and  the  MN  trajectory,  as  depicted 
in  Figure  0-2.  However,  as  the  number  of  unknown  variables  increases  and  they  become  more 
unpredictable,  the  SLAM  accuracy,  convergence,  and  hence  robustness  degrades. 
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Mapping  Tracking  _ | _ | 

SLAM  with  partial  SLAM 

knowledge 

Figure  0-2  Schematic  scenarios  of  a  SLAM  problem  in  a  broad  prospective 


Consider  a  MN  moving  from  left  side  of  Figure  0-2  where  the  MN  has  access  to  LOS 
signals  of  at  least  4  GPS  satellites  so  that  precise  trajectory  information  is  available. 
Opportunistic  wireless  localization  can  be  solved  for  the  estimation  of  the  stationary  ANs  (for 
example,  nearby  LTE  base  stations)  locations  or  the  mapping  problem.  As  the  MN  moves  into  a 
GPS-deprived  environment  such  as  an  urban  canyon  area  or  an  indoor  building  where  the  GPS 
signals  are  blocked,  it  has  no  access  to  the  precise  localization  information  but  has  prior 
knowledge  about  the  ANs  locations  from  the  previous  mapping  step.  The  OWLS  turns  to  the  MN 
tracking  problem  using  received  signals  of  opportunity  from  wireless  ANs  with  known  locations. 

As  the  MN  moves,  it  may  lose  sight  of  some  or  all  of  the  ANs  where  there  will  be 
partial/no  knowledge  about  their  source  locations,  as  depicted  in  the  right  side  of  Figure  0-2. 

This  is  the  worst  case  for  the  OWLS  since  it  must  not  only  map  the  source  locations,  but  also 
track  the  MN  trajectory.  Four  system  models  are  suggested  for  this  area,  as  addressed  in 
Chapter  3  and  Bayesian-based  solutions  are  proposed  as  the  final  stage  of  designing  the 
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tracking  algorithm,  as  addresses  in  Chapter  4.  In  first  system  model,  it  is  assumed  that  ANs  are 
synchronized  with  MN  and  also  all  ANs  are  synchronized  via  a  reference  clock  system  (like 
GPS  or  a  central  hub).  Moreover,  multipath  propagation  is  stationary  or,  if  not,  its  condition  is 
known  to  the  MN.  This  means  that  for  every  AN-MN  link,  the  propagation  is  either  LOS  or  NLOS 
in  the  whole  trajectory,  or  if  the  sight  condition  of  ANs  varies  during  trajectory  as  the  MN  moves, 
the  ANs’  LOS/NLOS  conditions  are  known  by  the  MN  via  prior  detection  processing,  and  NLOS 
measurements  are  discarded  from  the  observation  set. 

Second  system  model  considers  the  effect  of  MN  clock  drift  where  it  is  modeled  by  first- 
order  Markov  process  based  on  the  random  smooth  trajectory  assumption;  and  further,  in  third 
system  model,  it  is  assumed  that  MN  is  no  longer  synchronized  with  ANs  and  ANs  are  also  not 
synchronized  with  each  other.  This  assumption  introduces  a  range  offset  variable  for  every  MN- 
AN  link.  For  all  four  system  models,  both  EKF-based  and  PF-based  SLAM  solutions,  originally 
known  as  EKFSLAM  and  FastSLAM,  are  modified  for  newly-introduced  nuisance  parameters. 
Their  performances  are  compared  in  terms  of  consistency  and  efficiency  based  on  BFIM  and 
other  estimation  and  information  metrics  obtained  in  Chapter  3,  and  the  strengths  and 
weaknesses  of  each  method  are  discussed. 

In  the  forth  system  model,  the  ultimate  objective  of  this  thesis  is  addressed  that  is  to  develop  a 
SLAM-based  OWLS  for  the  scenario  where  not  only  the  MN  trajectory  and  AN  locations  are 
unknown,  but  also  where  LOS/NLOS  conditions  in  AN-MN  links  are  non-stationary  and 
unknown.  To  deal  with  different  propagation  conditions,  two  measurement  models  are 
introduced  for  LOS  and  NLOS  conditions,  where  the  transition  between  the  LOS  and  NLOS 
modeled  with  a  Markov  chain.  The  binary  state  of  LOS/NLOS  conditions  turns  the  OWLS  to  a 
jump  Markov  non-linear  system  (JMNLS)  [xiij.  Four  novel  SLAM-based  solutions  for  this 
nonlinear  JMNLS  are  suggested  combining  local  maximum  likelihood  (LML)  and  interacting 
multiple  model  estimators  with  EKF-based  and  PF-based  SLAM  solutions:  LML-EKFSLAM, 
IMM-EKFSLAM,  LML-FastSLAM  and  MM-FastSLAM.  IMM-EKFSLAM  is  derived  by  extension  of 
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the  EKF  for  joint  estimation  of  the  discrete  sight  states  and  continuous  variable  states  including 
the  MN  location,  AN  locations  and  range  error  using  interactive  multiple  model  (IMM)  estimator. 
In  MM-FastSLAM,  the  Rao-Blackewellized  (RB)  PF-based  SLAM  [xxxiii]  is  modified  based  on 
the  multiple  model  Particle  filter  (MMPF)  [xii],  where  the  posterior  PDF  samples  represent  both 
MN  location  dynamic  state  and  AN  sight  condition  states.  LML-FastSLAM  and  LML-EKFSLAM 
are  obtained  by  combining  modified  FastSLAM  or  EKFSLAM  proposed  for  system  model  3  with 
LML  criteria  for  detecting  ANs  sight  discrete  state.  Finally,  the  performance  of  each  method  is 
evaluated  with  the  related  BFIM  and  is  compared  with  previously  developed  methods  when  the 
NLOS  effect  and/or  range  error  are  ignored. 

Chapter  2  OWLS  Background  Studies  and  Challenges 

While  the  Global  positioning  system  (GPS)  is  widely  accessible,  its  position  accuracy  is 
compromised  in  dense  urban  areas  and  indoor  environments  due  to  low  signal  to  noise  ratio, 
LOS  blockage  and  narrowband  signaling  [xiii].  Such  factors  are  limiting  when  sub-meter 
accuracy  is  required.  Taking  advantage  of  opportunistic  signals  from  the  available  local  wireless 
networks,  such  as  locally  generated  WiFi,  WLAN  or  4G  LTE  wireless  signals  of  ample  power 
and  typically  larger  bandwidth  can  partially  ameliorate  these  issues.  The  wireless  signaling  is 
assumed  to  be  sourced  from  wireless  network  access  nodes  (ANs)  which  are  primarily  used  for 
mobile  data  and  voice  communications  [vii]. 

Similar  to  any  wireless  localization  method,  the  opportunistic  wireless  localization 
exploits  the  received  signal  observables  to  extract  the  MN  position  information.  In  addition  to 
unknown  MN  position  variables  in  typical  wireless  localization  methods,  the  tracking  algorithm 
adopted  for  the  opportunistic  wireless  localization  must  also  deal  with  the  estimation  of  unknown 
nuisance  parameters  of  signals  of  opportunity  from  opportunistic  sources  [xiv],  [iii].  Prior  to 
design  any  tracking  algorithm  for  an  OWLS,  it  is  required  to  present  a  fundamental  study  on 
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wireless  localization  methods,  opportunistic  wireless  localization  system  parameters,  and  its 
challenges. 

This  chapter  covers  the  required  background  studies  for  wireless  localization  based  on 
signal  of  opportunity  reception.  Primarily,  it  explains  the  wireless  positioning  methods  according 
to  signal  observables.  However,  the  wireless  signaling  is  assumed  to  be  sourced  from  local 
wireless  sources  such  as  4G  LTE  access  nodes  which  are  uncooperative  with  poorly  defined 
source  locations.  In  this  thesis,  the  opportunistic  sources  are  also  referred  as  uncooperative  or 
unregistered  network;  this  means  that  the  network  location  processing  or  assistance  is  not 
involved.  The  downlink  signaling  is  specifically  exploited,  as  the  MN  processing  is  desired  to  be 
capable  of  standalone  localization.  As  mentioned  earlier,  since  the  MN  receiver  is  not 
necessarily  registered  within  wireless  source  networks;  it  must  deal  with  the  detection  of  a 
multitude  of  signals  with  unknown  and  random  parameters  including  signal  characteristics,  such 
as  bandwidth,  the  carrier  frequency;  or  channel  properties,  such  as  shadowing,  multipath  and 
NLOS/LOS  conditions;  and  network  properties,  such  as  the  AN  location  and  clock 
synchronization.  In  this  regard,  we  also  discuss  the  opportunistic  wireless  localization 
challenges  due  to  signal,  channel  and  network  properties  in  more  detail.  Previous  studies  are 
also  presented  and  recommended  solutions  are  discussed.  Finally,  a  systematic  approach  to 
formulate  the  opportunistic  wireless  localization  problem  is  introduced  including  all  nuisance 
parameters  explained  in  former  sections. 

Signal  characteristics  used  for  positioning  include  the  angle  of  arrival  (AOA),  the  time  of 
arrival  (TOA)  or  the  time  difference  of  arrival  (TDOA),  the  signal  strength  (SS),  and  the  carrier 
phase  or  carrier  phase  difference  (also  known  as  the  phase  of  arrival  (POA/PDOA) ).  Traditional 
positioning  methods  utilize  triangulation  algorithms  that  exploit  the  geometric  properties  of 
triangles  to  estimate  the  MN  location  through  lateration  or  angulation. 

Angulation  estimates  the  location  of  a  MN  by  computing  the  angle  relative  to  multiple 
reference  points  via  directive  antennas  or  antenna  arrays.  The  accuracy  of  the  position 
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estimated  by  AOA  methods  is  very  sensitive  on  the  distance,  geometry,  and  multipath  in  the 
MN-AN  link  [xv],  In  angulation  or  AOA-based  approaches,  the  location  of  the  desired  target  can 
be  found  at  the  intersection  of  several  AOA  direction  measurements,  which  are  in  the  form  of 
ray  emanating  from  the  AN.  As  shown  in  Figure  0-3. a,  this  method  may  use  at  least  two 
reference  points  for  2D  position  information.  For  AOA  estimation,  no  synchronization  is  required 
between  transmitters.  However,  the  application  of  directional  antennas  or  antenna  arrays  in  a 
handheld  receiver  is  limited  due  to  the  complexity  and  size  of  antenna  arrays.  One  of  the 
disadvantages  of  AOA-based  positioning  methods  is  that  their  performance  is  very  sensitive  to 
accuracy  of  angle  measurements.  However,  high  accuracy  angle  measurements  in  indoor 
environments  is  limited  by  shadowing  and  multipath  reflections  arriving  from  misleading 
directions  that  may  be  as  strong  as  LOS  component  or  even  stronger,  or  by  the  directivity  of  the 
measuring  aperture. 


Figure  0-3  Geometric  locus  of  the  MN  location  estimate  in  wireless  positioning 

methods 
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Lateration  methods  (also  known  as  range-based  methods)  estimate  the  MN  position  by 
measuring  the  range  from  multiple  ANs.  Range-based  methods  utilize  the  SS,  TOA/TDOA, 
round-trip  time  of  flight  (RTOF)  or  POA/PDOA  measurements  of  the  received  signal.  TOA- 
based  algorithms  are  based  on  the  fact  that  the  distance  from  the  AN  to  the  MN  is  directly 
proportional  to  the  propagation  time.  As  shown  in  Figure  0-3. c,  the  position  is  computed  by  the 
intersection  point  of  the  circles  of  TOA  measurements.  One  of  the  disadvantages  of  these 
methods  is  that  all  of  the  ANs  and  the  MN  must  be  time  synchronized.  In  contrast,  TDOA-based 
algorithms  compute  the  relative  position  of  the  MN  by  examining  the  difference  in  TOAs  from 
various  ANs.  Here,  the  MN  location  lies  on  a  hyperboloid  with  a  constant  range  difference 
between  two  ANs,  as  depicted  in  Figure  0-3. d.  Unlike  TOA-based  methods,  they  only  requires 
ANs  to  be  synchronized. 

RTOF-based  methods,  similar  to  radar,  determine  the  time  of  flight  of  signal  traveling 
from  the  AN  to  the  MN  and  back,  as  illustrated  in  Figure  0-3. b.  For  these  methods,  a  more 
moderate  relative  synchronization  system  is  required  in  comparison  with  TOA  methods. 
However,  the  exact  delay/processing  time  caused  by  the  reference/responder  is  unknown  which 
cannot  be  ignored  in  short  range  indoor  applications. 

SS-based  or  signal  attenuation-based  methods  try  to  compute  the  signal  path  loss  due 
to  propagation  and  then  translate  it  into  the  range  estimate.  Due  to  the  dynamic  and 
unpredictable  characteristics  of  indoor  environment  propagation,  an  explicit  formulation  between 
distance  and  path  loss  does  not  always  hold.  For  this  reason,  SS-based  positioning  algorithms 
are  mostly  used  in  site-specific  applications  using  fingerprinting  and  scene  analysis  methods. 

Similar  to  time-based  methods,  carrier  phase-based  positioning  methods  are  considered 
as  range-based  positioning  methods.  However,  they  provide  better  position  accuracy  compared 
to  other  time-based  ranging  methods,  since  the  carrier  phase  can  be  measured  with  the 
precision  of  a  0.01-0.05  cycle.  Assuming  the  transmitting  signal  is  a  pure  sinusoid,  the  delay  can 
be  expressed  as  a  fraction  of  wavelength.  As  long  as  the  phase  is  in  interval  \0,2n)  ,  the  range 
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estimate  is  unambiguous;  otherwise,  it  is  necessary  to  deal  with  ambiguity.  GPS  receivers 
acquire  the  phase  lock  with  the  transmitter  signal,  measure  the  initial  partial  phase  difference 
between  the  received  and  transmitter  signals  (generated  at  the  receiver),  and  then  track 
changes  in  phase  difference  (i.e.  counting  full  cycles  and  keeping  track  of  the  partial  phase).  In 
an  ideal  LOS  condition,  the  integer  full  cycle  numbers  cannot  be  measured  and  have  to  be 
estimated;  however,  they  remain  constant  as  long  as  the  carrier  phase  tracking  loop  is  locked. 

Generally,  for  absolute  positioning  based  on  the  phase  measurement,  two  receivers  are 
required;  one  is  the  reference  receiver  at  the  stationary  known  location  and  the  other  is  the  MN 
receiver.  This  architecture  is  required  for  resolving  ambiguity  and  initialization.  After  initialization, 
the  MN  receiver  is  expected  to  track  the  phase  continuously.  The  integer  ambiguity  can  be 
resolved  using  the  least-squares  ambiguity  decorrelation  adjustment  (LAMBDA)  method,  carrier 
phase  measurements,  or  multiple  epochs  from  several  transmitters.  Since  the  integer  ambiguity 
solution  is  not  the  main  focus  of  this  thesis,  this  information  suffices  for  this  research  objective, 
and  the  interested  reader  is  referred  to  [xvij. 

The  mathematical  model  for  the  carrier  phase  measurement  in  units  of  cycle  in  an 
opportunistic  signal  reception  at  the  MN  receiver  is  described  as 

<p(t)  =  jr(t)  +  - StAN)  +  N  +  y>0  +  •>  °'2 

where  cp(t)  is  the  partial  carrier  phase  cycle  measured  by  the  receiver,  s  models  the 

measurement  noise,  and  (pQ  is  the  unknown  carrier  phase  of  the  AN  transmitter  at  the  time  of 

transmission.  The  A  and  f  are  the  carrier  wavelength  and  the  carrier  frequency,  respectively. 

r(t )  is  the  geometric  range  between  the  AN  and  MN  at  time  t.  The  MN  receiver  and  AN 

transmitter  clock  biases  are  represented  as  Stm  and  StAN ,  both  in  units  of  second.  N  is  the 

integer  ambiguity,  which  is  the  total  number  of  full  carrier  cycles  between  the  AN  and  the  MN.  It 
must  be  noted  that  since  the  range  measurement  from  the  carrier  phase  is  ambiguous,  it  is 
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referred  to  as  the  pseudo  range  in  many  contexts.  The  pseudo-range  measurement  is  defined 
as  a  function  of  the  true  range  measurement  which  usually  differs  by  an  unknown  offset. 

The  MN  is  not  necessarily  registered  in  the  network,  thus  the  ANs  do  not  provide  any 
synchronization  information  for  the  MN.  Consequently,  we  are  only  interested  in  carrier  phase 
changes  in  time  that  are  independent  of  the  initial  carrier  phase  at  the  AN  and  the  ambiguous 
integer  number  of  cycles.  Without  the  loss  of  generality,  we  assume  that  (p{t  =  0)  is  zero  that 
the  MN  tracking  is  initialized  at  the  coordinate  origin,  so  the  Eq.(O-l)  for  carrier  phase  variation, 
A cp(t) ,  as  a  function  of  time,  is  modified  by 

Ap(0  =  P(0-0  =  ^[K0-Kf  =  0  }  +  f.t(StMN  - dtAN )  +  sip  0-3 

where  St^  and  StAN  are  the  MN  and  the  AN  clock  rate  drifts  ,  respectively.  The  unknown  clock 
bias  term,  due  to  initial  non-synchronization,  is  eliminated,  since  it  is  considered  unchanged  in 
time  interval  [0,f] .  It  is  also  assumed  that,  in  the  time  interval  [0,f] ,  the  integer  number  of 
cycles  is  constant  and  there  is  no  cycle  slip  in  the  tracking  loop.  Figure  0-4  shows  general 
carrier  phase  relation  with  the  AN-MN  distance  in  LOS  propagation,  as  defined  in  Eq.  (0-2). 


AN 


MN 

Figure  0-4  Ranging  based  on  the  carrier  phase  measurement 
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Consider  a  small  portion  of  the  MN  trajectory  where  the  signal  is  collected  between  M 
sampled  positions,  p  =[p0,p15...,pM]r  ,  within  the  trajectory.  If  the  trajectory  length  is  within  the 


spatial  coherence  interval  of  the  channel  and  also  much  shorter  than  the  AN-MN  distance  r(t ) , 
the  plane  wave  propagation  over  the  M  sampled  locations  of  the  MN  trajectory  is  justified.  In  this 
regard,  the  signal  received  at  the  M  spatial  sampling  points,  p  =  [p0,p„...,pM_1]r,  is  denoted  as 


!(^P)  = 


’(^Po) 


0-4 


/(f.P«-i)_ 

Consider  the  coordinate  system  origin  at  the  initial  location  of  the  MN,  Eg.  (0-3)  in  terms 


of  time  delay,  is  written  as 

S(t-T0) 


i(Ap)  = 


0-5 


where  r  -  Pm 


.  is  the  unity  vector  of  direction  in  a  LOS  propagation  between  the 


AN  transmitter  and  the  MN  receiver,  defined  by 
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Assuming  that  5(f)  is  a  bandpass  signal,  we  have 

s(t)  =  V2  Rejs(/;)eM''},  m  =  -1 


0-7 


where  5(f)  is  the  lowpass  band  limited  complex  envelope  of  the  received  signal,  and  a>c 

is  the  carrier  frequency  in  radians.  Similarly,  the  signal  at  the  spatial  sampling  point  is 
modified  in  the  complex  envelop  representation  by 

0-8 

s (t, pm )  =  s(t-Tm)  =  J 2  Re {5 ( t - x m ) eJWc(,~Tm) } ,  m  =  0,...,M-1 
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If  it  can  be  assumed  that  incoming  signals  are  narrowband,  implying  that  the  reciprocal 
of  a  maximum  propagation  delay  across  M  MN  locations  is  much  greater  than  the  signal 
bandwidth, 


80 


A  T 
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where  Bs  is  the  bandwidth  of  the  complex  envelope  and  A7jrax  is  the  maximum  propagation 
delay  between  each  two  MN  locations  [xvii],  the  narrowband  assumption  justifies  as 


)□.?(?),  m  =  0,1,...,  M-l 


0-10 


Hence,  in  the  narrowband  case,  Eq.(0-7)  reduces  to 

s(7,pffl)  =  V2  Rej s(t)ei<°cte~Ja>cTm  | ,  m  = 
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In  a  narrowband  signal  model,  the  sensor  wise  propagation  delay  is  approximated  by  a 
phase  shift  as  a  function  of  the  propagation  direction,  signal  wavelength  ,  MN  position  (here,  the 
MN  displacement) 

0-12 

'Ik  j 

acrm  =Yk'W  Pm 

On  the  other  hand,  the  carrier  phase  measurement  in  units  of  cycle  for  positioning 
purpose  is  modified  as 

0-13 

1  T 

^(0  =  -k^  P (t)  +  (pclk+£y 

where  <pclk  is  clock  instability  due  to  free  running  clocks  in  receiver  (here,  the  MN)  and 
transmitter  (here,  the  ANs),  e  models  the  phase  error  due  to  both  variation  of  antenna  phase 
center  along  with  the  propagation  direction  and  phase  discriminator  noise. 
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Figure  0-5  Ranging  based  on  carrier  phase  measurement  in  a  plane  wave  propagation 

The  performance  of  OWLS  is  highly  dependent  on  both  the  quality  of  measurements  and 
unknown  characteristics  of  the  signal  which  interfere  as  nuisance  parameters  in  the  position 
estimation  problem.  In  this  regard,  it  is  important  to  explore  the  diverse  set  of  factors  that  can 
affect  localization  accuracy.  In  summary,  the  main  deficiency  factors  in  the  opportunistic 
wireless  localization  are  classified  as  unsynchronized  reception  from  an  uncooperative  network, 
non-stationary  wireless  propagation,  lack  of  knowledge  about  source  locations,  and  noise  and 
interference  from  unwanted  sources.  The  following  sections  discuss  in  more  detail  these  four 
typical  sources  of  error  in  opportunistic  wireless  localization. 

Indoor  localization  requires  precise  timing  measurement  in  at  least  nano-seconds  since 
the  accuracy  level  of  several  decimeter  is  desired.  For  this  reason,  clock  synchronization  is  an 
important  factor  in  the  ranging  and  positioning  performance.  The  carrier  phase  and  time-based 
positioning  methods  require  the  receiver  and  transmitter  nodes  to  be  equipped  with  a  stable 
oscillator  from  which  an  internal  clock  reference  is  derived  to  measure  the  true  time  with  high 
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accuracy.  However,  various  physical  effects  cause  oscillators  to  experience  frequency  and 
phase  drift,  leading  to  large  errors  in  phase  and  time  measurements,  and  therefore  in  the 
position  estimate.  In  a  disciplined  oscillator,  the  frequency  is  controlled  by  the  internal 
microprocessor  based  on  the  measurement  of  its  frequency  relative  to  the  received  GPS  or 
other  wireless  signals.  Any  interruption  in  signal  reception  may  degrade  the  clock  frequency. 

In  an  opportunistic  signal  reception,  while  the  MN  receiver  is  not  locked  to  the 
synchronization  signal  sent  from  the  AN,  an  improved  long-term  and  short-term  stabilized 
oscillator  is  of  great  importance.  The  clock  instability  characteristics  are  defined  by  two  terms, 
holdover  accuracy,  and  short-term  stability.  Holdover  accuracy  refers  to  the  time  error 
accumulation  during  a  free-running  clock  operation  in  a  long  term  observation  that  is  typically  1 
day.  For  example,  a  medium-stability  OCXO  accumulates  100  microseconds  of  time  error  in  24 
hours  after  the  reference  signal  (like  GPS  or  CDMA  signal)  is  lost.  Short-term  stability  refers  to 
the  frequency  deviation  relative  to  a  reference  frequency  standard,  measured  over  an 
observation  time  of  one  second  or  less.  When  measured  in  units  of  time,  it  is  known  as  jitter, 
and  when  measured  in  units  of  phase,  in  a  1  Hz  bandwidth  at  specific  frequencies,  is  defined  as 
phase  noise  [xviii].  It  implies  that  short-term  instability  characteristics  of  the  receiver/  transmitter 
clocks  mainly  affect  positioning  performance  rather  than  holdover  accuracy. 

The  local  time  of  a  clock  can  be  expressed  as  C(t )  as  a  function  of  real  time  t.  The 
oscillator  frequency  determines  the  rate  at  which  the  clock  runs.  In  a  perfect  clock,  the  rate 
denoted  as  dC{t ) / dt  would  be  equal  to  1 ,  or  C(t )  =  t .  However,  all  clocks  are  subject  to  clock 
drift  due  to  various  physical  sources  which  result  in  uncertainty  between  the  MN  and  the  rest  of 
the  ANs  network,  as  shown  in  Figure  0-6. 
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Figure  0-6  Illustration  of  ideal  time  and  local  clock  time  (short  term  instability) 

A  reasonable  assumption  that  will  be  made  here  is  that  the  ANs  network  is  tightly 
synchronized  perhaps  locked  to  a  central  synchronizing  clock  or  individually  to  a  GPS  source. 
This  is  a  rather  benign  assumption  to  make  as  each  AN  will  have  a  backhaul  connection  to  a 
communication  hub  of  sorts.  However,  the  actual  signal  delay  from  this  central  hub  to  ANs  is 
difficult  to  actually  control  and  maintain.  Regardless,  it  will  be  assumed  that  the  AN  clocks  are 
perfectly  stable  and  that  the  clock  instability  resides  with  the  MN  clock.  As  an  example,  the  LTE 
hub  is  likely  based  on  a  higher  quality  clock  than  the  consumer  grade  MN  clock  of  a  handset 
device. 

For  the  MN  free-running  clock  node  without  any  synchronization  scheme,  one  solution  is 
that  the  MN  receiver  can  adequately  model  the  behavior  of  its  own  clock.  The  time  reading  error 
of  a  clock  can  be  obtained  as  a  function  of  oscillator  quartz  frequency  given  by 


0-14 
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where  t0  and  C(t0 )  are  the  reference  time  epoch  and  clock  reading  at  the  reference  epoch; 
and  ft{t )  and  /0  are  the  frequency  of  the  oscillator  and  nominal  oscillator  frequency, 

respectively.  The  oscillator  frequency  shows  deterministic  linear  variations  as  well  as  random 
errors  [xix],  A  standard  model  for  the  frequency  of  an  oscillator  is  defined  as 


m=fb+f(c(t0)-t0)+f.(t) 


0-15 


where  both  fb  and  / ,  as  frequency  bias  and  frequency  drift,  model  the  linear  deterministic 


portion,  and  fr(t)  represents  the  unmodelled  random  frequency  errors. 

The  standard  approach  to  deal  with  clock  frequency  deviation  is  to  model  the 
deterministic  part  by  an  explicit  polynomial-like  function  and  model  the  random  part  by  its 
sample  variance  known  as  Allen  variance  [vii],  [xx],  [xxi],  [xxii],  Allen  variance  is  defined  as  the 

sample  variance  of  y(t )  =  ^  —  ,  given  by 

fo 


1  M  * 

*  m+T 

where  ym  =  -  J  y(t)dt 


0-16 


Allen  deviation  is  the  square  root  of  Allen  variance  that  is  regarded  as  an  indication  of 
the  overall  clock  stability  [xxiii],  [xxiv],  [xxv].  Figure  0-7  depicts  the  clock  stability  of  several  clock 
sources  in  term  of  Allen  deviation.  The  clock  stability  is  defined  as  a  function  of  observation  time 
of  a  particular  clock.  If  it  is  assumed  that,  at  the  start  of  observation  time,  the  clock  is 
synchronized  with  a  true  time  scale,  the  RMS  value  by  which  the  clock  has  deviated  after  a 
certain  time  interval  z  is  stated  by  r.crv(r)  . 
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Time  Domain  Stability  --  Oscillator  Options 


Tau,  seconds 


— • —  TCXO  MS-OCXO  — HS-OCXO  — m — US  OCXO 


— n—  HS-Rb 


Figure  0-7  Square-root  of  the  Allan  variance  of  typical  oscillators  (from 
http://www.endruntechnologies.com/pdf/OscOptions.pdf,  [xxvi]) 


As  an  example,  the  Allen  variance  of  10'9  for  a  TCXO  oscillator  implies  instability  of 
about  0.6  nsec  per  second  which  roughly  equates  to  20  cm  per  second  in  range  measurements. 
This  means  that  if  a  localization  scenario  consists  of  a  single  MN  and  AN,  and  the  location  is 
perfectly  known  at  t=0,  then  that  location  will  randomly  drift  as  a  random  walk  at  a  deviation  of 
about  20  cm  per  second.  However,  as  it  is  a  random  walk,  the  deviation  proportionally  increases 
to  the  square  root  of  the  time  lapse  in  seconds,  and  not  linearly  with  time.  From  the  Allen 
deviation,  it  is  possible  to  estimate  instability  in  terms  of  physical  distance  by  a  Markovian 
process. 

However,  mainstream  published  research  on  clock  instability  error  has  mostly  focused 
on  identifying  models  and  their  parameters  from  plots  of  Allen  variance  or  of  power  spectral 
densities;  whereas,  the  clock  stability  models  based  on  the  autocorrelation  function  of  phase 
noise  also  show  that  the  variation  of  clock  deviation  can  be  fitted  as  a  Markov  process,  [xxvii] 
presents  a  clock  deviation  modeling  based  on  a  second  order  random  process  where  the  model 
parameters  are  obtained  by  its  autocorrelation  function.  The  Markovian  process  of  clock  drift 
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helps  the  clock  deviation  parameters  to  be  modeled  along  with  the  dynamic  MN  trajectory  state- 
space  vector.  Consequently,  an  optimum  solution  for  the  OWLS  has  to  not  only  localize  the  MN 
position  as  it  undergoes  random  motion,  but  also  to  track  the  clock  deviation  to  mitigate  the 
localization  errors  due  to  clock  instability. 

The  main  source  of  error  in  any  wireless  localization  system  is  derived  from  multipath 
fading,  shadowing  and  LOS  blockage.  In  indoor  environments,  due  to  reflection,  detraction,  and 
scattering  of  radio  waves  by  structures  inside  a  building,  several  copies  of  transmitted  signal 
reach  the  receiver;  this  phenomenon  is  called  multipath.  In  narrowband  signaling,  different 
propagation  paths  are  added  constructively  (in-phase)  or  deconstructively  (out-of-phase)  within 
signal  propagation  time,  and  cannot  be  resolved.  The  multipath  delay  spreading  increases 
significantly  the  uncertainty  of  TOA  measurements  which  leads  to  large  bias  error  in  the  range 
estimate.  The  increased  bandwidth  partially  ameliorates  this  situation  because  the  channel  acts 
more  frequency  selective,  and  the  multipath  components  become  resolvable.  Therefore,  the 
leading  edge  of  received  signal,  corresponding  to  the  LOS  component  is  not  biased  with  the 
resolvable  multipath  returns.  As  an  example,  the  currently  deployable  LTE  sources  are  limited  to 
20  MHz  in  terms  of  bandwidth.  This  is  significantly  better  than  the  GPS  or  IS2000  signals  which 
are  typically  limited  to  about  1  MHz;  even  though,  this  is  not  sufficient  to  provide  accurate  TOA 
measurements  for  indoor  localization. 

LOS  propagation  may  also  be  obstructed  completely  or  partially  so  that  the  transmitted 
signal,  before  receiving  at  the  MN,  is  reflected  by  structures  or  travels  through  different 
obstacles  like  walls  in  a  building.  Both  of  these  phenomena  lead  into  a  positive  random  excess 
delay  in  time  measurements  that  render  traditional  methods  of  localization  virtually  ineffective.  In 
other  words,  combating  the  effects  of  multipath  is  the  most  significant  challenge  in  providing 
adequately  accurate  position  estimates,  and  hence  new  methods  are  required. 

In  this  regard,  many  research  works  have  focused  on  LOS/NLOS  identification  as  prior 
processing  of  any  localization  algorithm.  The  proposed  methods  to  mitigate  the  NLOS  effect  are 
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categorized  in  two  groups.  In  the  first  group,  the  NLOS/LOS  conditions  are  detected  in  a  prior 
processing  to  localization.  Whenever  the  LOS/NLOS  conditions  are  identified,  NLOS 
measurements  can  simply  be  ignored  or  the  measurement  model  is  adopted  based  on  the 
range  error  to  improve  accuracy.  In  a  mixed  scenario,  the  carrier  phase  of  the  AN  transmission 
will  be  used  when  LOS  conditions  are  evident;  in  NLOS  situations,  the  code  phase  is  used  with 
multipath  delay  averaging. 

In  the  second  group,  sight  condition  and/or  its  related  parameters  are  jointly  estimated 
and  tracked  along  with  location  states.  In  a  Jump-Markov  system,  the  algorithm  must  solve  for 
both  discrete  (here,  binary  state  regarding  NLOS/LOS  conditions),  and  continuous  (here, 
position  variables)  states  at  the  same  time.  Jump  Markov  or  switching  Markov  model  is  also 
extended  for  time-based  systems  where  the  range  measurements  in  NLOS  condition  is 
modeled  by  jump  in  variance  and  mean  of  measurement  noise  [xlix-li],  [Ivj.  In  NLOS  situations, 
the  MN  motion  can  be  used  to  accumulate  a  number  of  TOA  measurements  from  the  AN,  and 
hence  averaging  out  the  effects  of  multipath.  This  is  essentially  handled  by  the  Bayesian 
processing  which  has  enhanced  indoor  positioning  performance  in  NLOS  situations. 

While  the  MN  receiver  must  deal  with  detection  and  tracking  of  the  signal  of  opportunity 
from  a  local  wireless  network,  it  has  to  combat  the  interference  of  unwanted  signals  from  other 
coexisting  devices  using  the  same  frequency  band.  The  OWLS  is  expected  to  operate  in  the 
presence  of  both  interferences  from  narrowband  interferes  such  as  Bluetooth  devices  and 
cordless  phones  and  wideband  interferes  such  as  microwave  ovens  and  WiFi  devices. 
Beamforming  methods  [xxviii],  [xxix]  and  nonlinear  filtering  are  considered  as  effective 
approaches  to  mitigate  the  effect  of  interference  on  localization  [xxx],  [xxxi],  [xxxiij.  However, 
the  effect  of  interference  on  OWLS  and  its  mitigation  techniques  is  out  of  the  scope  of  this  thesis 
and  in  the  proposed  OWLS  modeling,  the  interference  is  simply  treated  as  a  part  of  system 
noise  into  a  spread  spectrum  MN  device. 
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The  overall  objective  of  a  tracking  problem  is  to  generate  the  belief  map  of  the  MN 
location  as  the  trajectory  unfolds.  The  belief  map  will  provide  a  posterior  estimate  of  the  MN 
location  that  is  iterated  at  each  time  step  based  on  all  available  data  up  to  the  current  time.  The 
Bayesian  filtering  methods  such  as  PF  and  KF/EKF  provide  an  incremental  estimate  of  the 
posterior  PDF  based  on  the  accumulation  of  the  MN  trajectory  estimate  and  the  measurements 
made  up  to  current  time  with  knowledge  of  AN  locations. 

However,  in  an  opportunistic  signal  reception,  it  is  not  guaranteed  that  the  location  of 
ANs  are  known  to  the  MN.  This  gives  the  impetus  to  employ  the  SLAM  techniques  for  the 
OWLS  as  a  systematic  approach  to  take  into  account  all  available  constraints  regarding  the 
large  set  of  unknowns  such  as  the  MN  trajectory  and  AN  locations.  The  following  sections 
provide  a  short  history  of  the  SLAM  problem  as  well  as  its  classifications  and  proposed 
solutions.  Later,  observability  of  the  SLAM-based  OWLS  is  discussed. 


Origin  of  SLAM 

Simultaneous  Localization  and  Mapping,  which  is  also  known  as  Concurrent  Mapping 
and  Localization  (CML),  addresses  problems  where  the  MN  does  not  have  access  to  the  map  of 
the  environment,  nor  does  it  know  its  own  location.  SLAM  will  search  for  the  possible  tracking 
solution  for  a  MN,  if  placed  at  an  unknown  location  in  an  unknown  environment,  as  it 
consistently  builds  the  map  of  environment;  while  simultaneously,  estimates  its  location  within 
the  map. 

The  roots  of  SLAM  date  back  to  1986  IEEE  Robotics  and  Automation  Conference  in  San 
Francisco  when  probabilistic  mapping  as  a  fundamental  robotic  problem  is  addressed  and 
discussed  by  Peter  Cheeseman,  Jim  Crowley,  and  Hugh  Durrant-Whyte,  Raja  Chatila,  Oliver 
Faugeras,  Randal  Smith  [xxxiiij.  These  discussions  were  mainly  focused  on  application  of 
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estimation-theoretic  methods  in  Artificial  Intelligence  and  Robotics.  Subsequently,  many 
research  papers  have  come  out  on  related  problems  such  as  Smith  and  Cheesman  and 
Durrant-Whyte  [xxxiv]  Ayache  and  Faugeras  [xxxv]  ,  Crowley  [xxxvi]  and  Chatila  and  Laumond 
[xxxvii].  The  key  result  of  these  research  studies  proved  that  there  is  a  high  degree  of 
correlation  between  the  estimates  of  different  landmarks  locations  in  a  map,  which  increase  with 
future  observations  [xxxviii]. 

Later  in  1990s,  a  key  research  on  convergence  issues  was  developed  by  Csorba  [xxxix] 
and,  then  specifically,  on  Kaman-filter-based  SLAM  methods  and  the  probabilistic  localization 
and  mapping  methods  by  Thrun  [xl].  The  conceptual  breakthrough  in  SLAM  research  was  found 
that  the  combined  mapping  and  localization  problem,  once  formulated  as  a  single  estimation 
problem,  is  actually  convergent.  Most  importantly,  it  was  recognized  that  the  correlations 
between  landmarks  were  actually  the  critical  part  of  the  problem.  In  other  words,  the  higher 
these  correlations,  the  better  the  solution. 

Today,  the  SLAM  is  introduced  for  many  applications  in  a  number  of  different  domains 
from  indoor  robots  to  outdoor,  underwater,  and  airborne  systems.  It  has  been  formulated  and 
solved  as  a  theoretical  problem  in  many  different  forms.  The  SLAM  problem  is  a  continuous  and 
discrete  estimation  problem.  Generally,  the  continuous  estimation  part  pertains  to  both  the 
location  of  the  object  in  the  map  and  the  mobile  node’s  own  location,  and  the  discrete  nature  is 
related  with  correspondences.  When  an  object  is  detected,  a  SLAM  algorithm  must  distinguish 
the  relation  of  this  object  to  the  previously  detected  objects. 

Consider  a  MN,  such  as  a  vehicle  or  a  person  carrying  a  handset  device  moving  through 
an  unknown  environment  and  taking  observations  from  opportunistic  sources.  Lets  define  the 
following  variables  in  regards  to  the  OWLS  formulation: 

p(  :  the  state  vector  describing  the  location  of  the  MN  at  time  t . 
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u? :  the  control  vector,  if  available,  applied  at  time  t- 1  to  drive  the  MN  location  to  p,  at 
time  /;  or  observation  vector  that  can  model  the  MN  motion  state  process  such  as  computer 
vision  outputs. 

m; :  the  state  vector  describing  the  location  of  the  /th  ANs  whose  true  location  is 
assumed  to  be  time  invariant. 

m  :  the  stacked  vector  of  AN  locations,  m  =  {mi,  (  =  1,..,^}  where  Nan  is  the  number 

of  ANs. 


map :  the  stacked  vector  of  ANs  location  whose  locations  are  known  to  the  MN,  referred 
by  anchor  points  (APs),  m  =  {m;,  i  =  \,...,NAP) ,  where  N AP  is  the  number  of  APs 

m/;; :  the  stacked  vector  of  ANs  location  whose  location  are  unknown  to  the  MN,  referred 


by  feature  points  (FPs),  m  fp  =  {m;,  /  =  N AP  +1  ,...,NAP  +  NFP} ,  where  NFP  is  the  number  of  FPs. 
z;  t :  the  observation  vector  received  by  the  MN  from  /th  AN  at  the  time  t . 
zt  D  |z it,i  =  1,...,/^} ;  the  stacked  vector  of  observations  received  by  the  MN  from  all 
ANs  at  the  time  t . 

In  addition,  the  following  sets  are  also  defined: 

p1:(  □  {pj  p2  ...  P,}-{Pi:,_i  p(} ;  the  history  of  the  MN  locations. 


U17  ={U1  U2 

ZM:f  B  {Zi,l  Z/,2 


uj  =  {uH_j  uj  ;  the  history  of  update  control  inputs. 
z; , }  =  |z ■  lt_j  z ;  the  history  of  observation  received  by  the  MN 


from  /'-th  AN. 

It  must  be  noted  that  Nan  =  NAP  +  NFP .  The  optimum  solution  for  a  SLAM-based  OWLS 

must  develop  a  systematic  way  of  accounting  for  the  all  observations  to  compute  the  posterior 
belief  map  of  state  variables  including  the  MN  location,  FP  locations,  and  nuisance  parameters 
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(if  available).  According  to  posterior  PDF  computation  process,  the  SLAM  techniques  are 
classified  into  two  categories:  Online  SLAM  and  full  SLAM.  The  online  SLAM  problem  includes 
estimating  the  posterior  PDF  over  the  current  MN  location,  p, ,  along  with  the  unknown  source 

location  mfp  (and  other  nuisance  parameters,  if  available)  represented  as  p( p,,m/;i  |  zw,uw) 

[xxxiii].  Online  SLAM  algorithms  are  incremental,  and  they  discard  past  measurements  and 
update  controls  that  have  been  processed.  The  second  category  is  known  as  the  full  SLAM.  The 
full  SLAM  involves  calculating  the  posterior  PDF  over  the  entire  path  plt  along  with  the  FP 

locations  m^,  i?(p1:f, |z1:„u1:f). 

Figure  0-8  illustrates  the  graphical  schemes  for  states  dependency  of  a  general  SLAM 
problem  and  also  the  difference  between  the  full  SLAM  and  online  SLAM.  The  arrows  define  the 
dependency  between  processes  and  dark  grey  boxes  show  the  unknown  states  which  must  be 
estimated  at  time  step  t .  The  arrows  connectivity  between  pA ’s  shows  that  given  the  control 

data  at  time  t ,  u,  (if  available),  the  MN  location  at  time  step  t  is  only  dependent  to  MN  location 
at  previous  time  t- 1 .  In  other  words,  the  MN  location  variables  pA  are  dependent  is  a  specific 
way  that  the  past  and  future  are  dependent  only  through  the  present.  Moreover,  it  is  seen  that 
even  though,  state-dependent  measurement  process  zt  evolves  independently  of  each  other 

and  are  only  dependent  on  random  process  p, ,  since  the  FP  map  m/p  is  unknown,  estimate  of 
p,  and  mfp  are  not  independent.  Consequently,  the  joint  posterior  PDF  p,  and  mfp  cannot  be 

separated.  The  online  SLAM  and  full  SLAM  are  fundamentally  different  algorithms  as  the  set  of 
variables  to  be  estimated  is  different  in  each  algorithm  as  illustrated  in  Figure  0-8.  However, 
online  SLAM  can  be  modified  as  the  integration  of  the  full  SLAM  over  past  locations. 
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Figure  0-8  Graphical  schemes  of  the  a)  online  SLAM  in  comparison  with  the  b)  full 

SLAM 

Estimation  error  correlation  effect  on  SLAM 

One  of  the  important  aspects  of  SLAM  is  to  understand  the  correlation  between  estimate 
of  the  MN  trajectory  and  the  estimate(s)  of  unknown  parameters  of  environment  map.  For 
OWLS  with  no  nuisance  parameters  except  FP  locations,  Figure  0-9  shows  the  estimates  of  the 
MN  trajectory  and  the  FP  locations  compared  to  the  actual  MN  trajectory  and  FP  locations.  It 
illustrates  that  the  estimation  error  for  the  FP  locations  is  highly  correlated  with  the  MN  locations 
estimation  error.  This  is  because  the  main  source  of  the  estimation  error  originates  from  the 
uncertainty  about  the  MN  position  where  the  observations  are  made  [xlj.  This  common  source  of 
error  makes  estimation  error  of  the  feature  points  highly  correlated  as  it  shapes  the  joint 
probability  density  of  feature  points  to  be  highly  peaked,  even  if  the  marginal  density  of  each 
individual  FP’s  PDF  has  a  skewed  shape. 
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Figure  0-9  Correlation  between  the  estimation  error  of  the  MN  trajectoryand  FP  locations 
compared  with  the  true  MN  trajectory  and  the  FP  locations 


The  interesting  point  is  that  the  correlation  between  the  estimates  of  FP  locations 
increases  in  time  and  makes  the  peak  of  the  joint  PDF  sharper  as  more  observations  are 
received.  As  discussed  earlier,  this  correlation  effect  originates  from  the  constraints  defined  by 
disparate  sources  of  information  in  OWLS  such  as  MN  motion  process  and  observables  from 
unknown  or  known  source  locations.  As  already  illustrated  in  Figure  0-1 ,  the  effect  of  correlation 
between  the  MN  trajectory  estimates  and  FP  locations  cast  the  SLAM-based  OWLS  as  a 
network  of  springs  where  minimum  energy  solution  of  this  spring  system  solves  for  the  MN 
trajectory  and  FP  locations  (and  other  nuisance  parameters,  if  available).  An  observation  at 
each  time  step  is  analogous  with  a  displacement  in  the  spring  network,  and  its  correlation  effect 
on  its  neighbors  is  proportional  to  their  distances  to  other  feature  points,  i.e.,  the  nearer  they 
are,  the  greater  is  the  effect,  as  shown  in  Figure  0-10.  As  the  MN  moves,  and  receives  more 
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observations  from  feature  points,  the  correlation  between  the  estimates  of  feature  points 
increases  and  the  spring  becomes  stiffen  The  thicker  link  in  Figure  0-10  denotes  more  stiffness 
and  thus  higher  correlation. 

FP 

FP 


Figure  0-10  Analogy  of  the  estimation  error  correlation  in  a  SLAM  problem  with  a  spring 

network 

In  theoretical  limit,  the  rigid  or  accurate  relative  map  of  ANs’  location  is  obtained  and  the 
relative  MN  location  accuracy  becomes  equal  with  the  MN  location  accuracy  achievable  with  a 
given  map.  The  minimum  energy  solution  for  this  spring  network  is  the  function  of  origin  of 
coordinate  systems  and  its  unit  vector.  If  there  is  no  constraint  to  an  absolute  point  or  coordinate 
system,  as  illustrated  in  Figure  0-1 1  .a,  there  is  only  constraints  to  the  origin  of  the  coordinate 
system  but  not  the  unit  vectors,  as  illustrated  in  Figure  0-1 1  .b,  there  are  a  family  of  lowest 
energy  solutions  which  gives  rise  to  the  ambiguity  relative  to  absolute  space.  The  correlation  is 
a  result  of  the  minimum  energy  solution  mapping  to  such  a  family  of  solutions  that  are 
ambiguous  as  far  as  an  affine  transformation.  That  is  the  family  of  minimum  energy  solutions  is 
a  set  of  solutions,  wherein  an  arbitrary  affine  transformation  exists  between  its  members. 

In  other  words,  the  location  accuracy  of  the  MN,  relative  to  the  map,  is  bounded  only  by 
the  quality  of  the  map  and  the  relating  measurement  sensors,  whereas  in  the  theoretical  limit,  it 
can  be  claimed  that  the  MN  relative  location  accuracy  becomes  equal  to  the  localization 
accuracy  achievable  with  a  given  map. 
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Consider  a  MN  that  moves  with  an  unknown  trajectory  and  receives  N  measurements 
at  each  time  step  from  N  FPs,  recalling  that  a  FP  is  an  AN  with  an  unknown  location  to  the  MN. 
For  a  minimal  2D  OWLS  where  there  is  no  nuisance  parameters  except  FP  locations,  there  is  at 
least  2  +  2 N  unknown  variables,  which  is  less  than  the  number  of  the  known  constraints,  N . 
Under  no  assumption  on  ANs  or  MN  trajectory,  after  K  time  steps,  the  number  of  unknown 
variables  increases  to  2 K  +  2 K x N  ,  while  there  are  only  KxN  known  equations,  leading  to 
an  unobservable  SLAM  problem.  The  first  assumption  to  make  a  SLAM  problem  observable  is 
stationary  ANs,  meaning  that  there  will  be  enough  time  steps  like  K  ,  where  2 K  +  2N  <  NK  or 
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K  >  — — — .  The  effect  of  stationary  ANs  assumption  for  a  two-time  step  SLAM  problem  with 


two  ANs  is  depicted  in  Figure  0-12.  It  can  be  seen  that  under  stationary  ANs  assumption,  four 
unknown  AN  location  states  are  merged  to  two,  decreasing  the  number  of  unknowns  to  8  for  2D 
localization  problem. 


AN  02 


Figure  0-12  Stationary  ANs  in  a  SLAM  problem 

If  the  MN  moves  in  an  agile  way  or  the  distance  between  measurements  is  so  far  that 
there  is  no  correlation  between  the  MN  positions,  no  information  related  to  the  previous  MN 
location  estimate  is  transformed  for  the  next  time  step  estimation.  The  problem  becomes  worse 
when  the  SLAM  problem  also  includes  an  estimation  of  random  nuisance  variables  such  as 
range  errors.  After  K  time  steps,  the  number  of  unknowns  is  increased  to  2 K  +  KN  +  2N , 
since,  at  each  time  step,  a  new  bias  variable  is  added  for  every  AN-MN  link,  while  there  are  only 
KxN  measurements  available.  In  this  case,  the  SLAM  problem  becomes  unobservable  unless 
we  can  assume  a  smooth  trajectory  for  the  MN  motion.  The  smooth  trajectory  assumption 
proposes  correlation  between  subsequent  dynamic  states  such  as  the  MN  position  and  range 
biases,  i.e.,  the  previous  or  future  states  estimate  can  add  information  for  the  current  state 
estimate,  in  addition  to  the  current  measurement.  This  information  is  conveyed  by  the  update 
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information  from  previous  states  in  the  online  SLAM,  or  from  previous  and  future  states  in  the 
full  SLAM. 


Bayesian  based  SLAM 

In  probabilistic  form,  the  optimum  estimate  for  an  online  SLAM  problem  is  obtained  by 
evaluating  the  posterior  PDF  p(qt  |  z1:,,u1:„q0) ,  where  q,  =[qd>f  q  jr  is  the  vector  of  unknown 

variables  including  dynamic  variables  qrf  t ,  such  as  the  MN  location  at  time  t ,  p, ,  and  random 
nuisance  parameters,  such  as  range  error,  and  stationary  variables  qv  ,  such  as  the  FP 
location,  mf  .  Joint  posterior  densities  of  FP  locations,  MN  location  and  other  nuisance 
parameters  embodies  all  the  statistical  information  that  can  be  extracted  from  observables  zu , 
control  data  ult  ,  and  prior  knowledge  q0  about  the  state  q, .  A  recursive  solution  for  the  SLAM 


problem  is  obtained  using  Bayes’  rule,  given  by 

P\Ht  I  L\\f>  Ul:*’4o/ 


0-17 


p(zt  I  Z  1:1-1’  Ul:r) 

Based  on  the  first  Markov  assumption,  as  shown  in  Figure  0-6,  given  q, ,  zt  is 
dependent  on  previous  measurements  and  control  data;  hence  Eq.  0-16  can  be  re-written  as 

n(n  „  jp(z,|q<)p(qjz1:r_1,u1„q0)  0-18 

p(zt\zl,-l’»l,) 

p{ zt  |  q,)  is  the  likelihood  PDF  that  describes  the  probability  of  making  an  observation 

when  FP  locations,  current  MN  location  and  nuisance  parameters  are  known.  Based  on  smooth 
random  trajectory  assumption,  the  dynamic  state  transition  is  considered  as  a  first-order  hidden 
Markov  process,  thus  the  prior  PDF p( q,  |  zl:(_pub,q0)  can  be  calculated  by 

f  0-19 

p(  q,  I  zi,-Pui^q0)  =  J  p^d,t  I  qrf,,-PuXq^-pq,  I  zi:f-Pup-pqo)^q^-i 
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Eq.  (0-17)  and  (0-18)  provide  a  recursive  procedure  for  calculating  the  joint  posterior 
PDF  based  on  history  of  observations  and  control  data  up  to  time  step  t .  Since  the  denominator 

p{ z,  |  z,M,uw)  in  Eq.  (0-17)  is  independent  of  q, ,  it  can  be  replaced  by  a  normalizing  constant, 
// ,  and  posterior  PDF  can  be  only  modified  in  terms  of  prior  and  likelihood  PDF: 

i7(qJzi^Ui„,qo)  =  ai7(zJqr)^(qJZi:(_i,Ui:„q0)  0-20 

^iKzJqMqJzi:r-i>ui:,>qo) 

The  optimal  Bayesian  estimate  of  the  state  q,  can  be  obtained  from  p(pnmf>  |zH)  by 

using  the  maximum  a  posterior  (MAP)  or  the  minimum  mean-square  error  (MMSE)  criterion. 
However,  in  the  general  case  of  nonlinear  dynamic  systems,  as  the  SLAM-based  OWLS 
problem,  the  posterior  PDF  cannot  be  determined  analytically  from  Eq.  (0-17)  and  (0-18). 
Particularly,  it  must  be  noted  that  the  dynamic  portion  of  the  state  vector  (including  the  MN 
location)  and  the  stationary  portion  of  the  state  vector  (including  FP  locations)  are  correlated 
through  measurements,  and  it  worth  nothing  to  compute  the  conditional  density  of  the  stationary 

states  p{(\s  |zw,ult,q0) ,  or  the  MN  location  p(qdl  |  zw,uw,q0) ,  separately,  since  both  are 

conditioned  on  the  deterministic  knowledge  of  other. 

Suboptimal  solution  to  the  SLAM  problem  requires  appropriate  presentations  of  the 
measurement  and  transition  processes  based  on  reasonable  assumptions  that  allow  the 
evaluation  of  the  posterior  PDF  with  efficient  computation.  The  EKF-SLAM  and  FastSLAM,  as 
suboptimal  Bayesian  filters  for  SLAM  problems,  are  proposed,  whereby  the  former 
approximates  the  posterior  PDF  with  a  Gaussian  PDF,  and  the  latter  describes  the  MN  states  as 
set  of  samples  (known  as  particles)  with  a  general  non-Gaussian  PDF. 


Extended  Kalman  Based  SLAM 

The  most  common  solution  to  a  SLAM  problem  in  the  form  of  a  state-space  model  with 
additive  noise  is  known  as  the  extended  Kalman  SLAM  (EKFSLAM).  The  first  element  and  the 
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basis  of  EKFSLAM  is  to  describe  the  evolution  of  the  dynamic  states  as  a  first-order  Markov 
process,  in  form  of 

0-21 

q^  =  /(q^-Pu,)+v. 

where  /(*)  is  a  known  deterministic,  possibly  a  nonlinear  function  of  previous  states  (\d  t  ]  , 
and  v(  □  p{ (\dt  |  qdM,u,)  is  additive,  zero  mean  ,and  uncorrelated  Gaussian  process  with 
covariance  Qt . 

The  second  element  in  the  EKF-SLAM  is  the  observation  model,  which  defines  the 
relation  between  observables  and  unknown  state  variables,  given  by 

z,  =//(q,)  +  w,  0-22 

=  Mq(/,„q.s)+w, 

where  h(*)  is  a  known  deterministic,  possibly  nonlinear  observation  function,  and 

wt  0  p( zt  |  qrf/,qi)  is  an  additive,  zero  mean,  and  uncorrelated  Gaussian  observation  noise 

with  covariance  matrix  Rt. 

The  standard  EKF-SLAM  algorithm  computes  the  estimation  based  on  the  mean  of  the 
conditional  posterior  PDF  />(q^,qs  |  z1:(,u1:<  j  as 


_ 
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q. 

and  covariance  matrix  of  estimation  is  defined  by  the  posterior  PDF  covariance  matrix  : 
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The  first  stage  of  EKF-SLAM  is  the  time-update  of  dynamic  states  or  prediction  stage. 


<1.,,  =  /(q, 

=V/T„j,MMv/r  +  a 

where  V/is  the  Jacobian  of  /(*)  evaluated  at  the  previous  estimate  . 
The  second  stage  is  the  observation  -update  stage: 

+  K,  [z,-ft(q,,.,|„q,.,)] 

S,„  =  I„_,  +K,D,K,r 
where 


q  d,t\t 
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0-27 


D  ,=Hl.itHHT+Rt 
K,  =  Ht D  1 

H  is  the  Jacobian  of  h  ,  evaluated  at  the  predicted  state  =  [qrf w_,]7  .  In  a 

EKF-SLAM,  if  the  determinant  of  covariance  matrix  monotonically  converges  to  zero,  it  is 
defined  that  EKFSLAM  is  convergent  in  stationary  states.  However,  the  convergence  is  only 
guaranteed  in  a  jointly  linear  Gaussian  system.  It  must  be  noted  that  standard  deviation  of  an 
individual  FP  estimate  converges  toward  a  lower  bound  imposed  by  initial  uncertainties  of  the 
MN  location  and  observations. 

Similar  to  EKF,  EKFSLAM  employs  a  linearized  approximation  of  the  nonlinear  motion 
and  observation  models;  this  inherits  many  caveats.  Nonlinearity  can  be  a  significant  problem  in 
EKFSLAM,  which  leads  to  inevitable  and  sometimes  drastic  inconsistency  in  solutions  [xxxiii], 
[xliv], 

FastSLAM 

Online  SLAM  methods  based  on  EKF  can  handle  the  large  number  of  variables  as  the 
update  of  the  posterior  PDF  only  involves  the  mean  vector  and  the  covariance  matrix.  However, 
some  of  the  state  variables,  notably  the  MN  location  variables  are  poorly  approximated  by  a 
joint  Gaussian  PDF  where  the  measurement  function  represents  a  highly  nonlinear  relation  with 
the  state  vectors.  Moreover,  non-Gaussian  prior  belief  maps  cannot  be  accurately  incorporated. 
PFs  or  sequential  Monte  Carlo  (SMC)  methods  is  suggested  as  an  alternative  implementation  of 
Bayesian  filters,  where  EKF  fails,  since  it  allows  for  a  general  PDF  of  the  state  variables 
conditioned  on  non-Gaussian  nonlinear  observables.  In  a  SLAM-based  OWLS,  the  PF  needs  to 
represent  all  state  space  variables  such  as  the  MN  trajectory  and  the  feature  point  locations,  by 
particles.  The  high  dimensional  state-space  of  the  SLAM  problem  makes  the  complete 
application  of  particles  infeasible  due  to  the  high  computational  complexity. 
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An  effective  solution  for  this  problem  is  the  Rao-Blackwellized  particle  filters  (RBPF). 
Rao-Blackwellization  which  is  also  known  as  variance  reduction  methods  are  based  on 
partitioning  the  state  vector  so  that  one  component  of  the  partition  is  conditionally  linear 
Gaussian  state-space  model  [xlij.  For  this  component,  analytical  solution  can  be  found  by 
exploiting  the  KF  while  the  PF  is  only  used  for  non-linear  non-Gaussian  portion  of  state-vector 
[xii], 

A  key  feature  to  apply  Rao-Blackwellization  theorem  in  a  SLAM  problem  is  that  the  FP 
constraints  become  mutually  independent  when  the  MN  trajectory  variables  are  given.  More 
precisely,  when  all  of  the  mutually  shared  variables  are  given,  then  the  remaining  FP  variables 
are  independent.  In  this  way,  an  EKF  can  be  applied  to  each  of  the  FP  variable  sets.  For  better 
understanding  of  FastSLAM  implication,  we  limit  dynamic  states  to  MN  location  variables  and 

stationary  states  to  FPs  location  variables,  where  q,  =  [p,  m/;Jr  ,  however  the  partitioning 
does  not  preclude  the  possible  augmentation  of  additional  variables. 

In  other  words,  the  i  -th  FP  location  vector,  m(  =  [mi  v  mi  v]r ,  is  independent  of  other 

FPs  location  variables,  when  the  mutually  shared  variables  including  the  MN  position  are  given. 
Hence,  the  joint  posterior  PDF  conditioned  on  the  set  of  measurements  and  deterministic 
updates  can  be  factorized  as  follows: 

Nam  0-28 

P\Pl:fmjp  I  Z1:,5U1:/)  =  /?(Pl:,  |  ZW,U1;(  )  Yl  H  lZl:,>Ul:,) 

‘=Nap+1 

The  product  term  indicates  the  independence  of  FP  state  variables,  when  conditioned  on 
the  MN  position  state  variables,  which  are  the  part  of  the  particles.  There  are  a  total  of  Np 

particles,  with  each  particle  consisting  of  the  complete  MN  trajectory  and  a  potential 
arrangement  of  the  FPs.  For  m  -th  particle,  we  have  the  trajectory  variables  of  p1”1  =  {x,yf™]  , 

and  the  Gaussian  FP  parameters  of  M[nIPnL]  ■  For  each  particle,  a  separate 

EKF  is  applied  for  each  FP  location  state.  Therefore,  multiple  possible  MN  trajectory  solutions 
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are  tracked  instead  of  only  the  most  probable  one.  The  FastSLAM  algorithm  stages  can  be 
described  as 


Retrieval  -  obtain  the  MN  location  of  the  previous  time  step  from  pj"11, 


Prediction  -  draw  new  MN  location  variables  based  on  the  proposal  distribution 


p["]Q  ^(pJp 


Measurement  update  -  update  FP  location  estimate  conditioned  on  the  MN  location. 

For  each  particle,  an  EKF  update  is  performed  over  the  observed  ANs  as  a  simple  mapping 

operation  under  known  MN  location  assumption.  The  outputs  of  this  stage  are  for  all 

particles  and  FPs  locations. 

PF  weight  update  -  calculate  the  importance  weight  wj"'] ,  which  is  based  on  the 
conditional  probability  of  p{zt  |  )  where  z,  is  the  set  of  measurements  from  all  ANs. 

(  \  N 

Resampling  -  sample  with  the  replacement  of  Np  particles,  from  the  set  {p[0"^]  |  r  , 
including  their  associated  maps,  considering  each  particle  has  probability  proportional  to  wj"1 . 

Given  Np  particles  and  Nfp  conditionally  indopondont  FPs,  thoro  will  bo  NpNpp  EKFs, 
each  with  a  small  cluster  of  variables  pertaining  to  the  specific  FP.  Consider  the  details  for  the 
EKF  of  the  m  -th  particle  and  the  i  -th  FP.  The  output  of  this  stage  is  .  defined  by 


0-29 


where,  for  i  -th  FP, 


0-30 


vl>]  _  y[w] 
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and 
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0-31 


where  R.t  is  the  covariance  matrix  of  the  measurement  set  zit  from  the  ith  FP. 

Once  the  EKF  is  finished  updating,  the  particle  weights  can  be  determined.  The  set  of 
each  FP  measurement  is  assumed  to  be  independent  and  the  likelihood  probability  is 
determined.  The  likelihoods  of  the  set  of  measurements  can  be  multiplied,  forming  the  overall 
weight  of  the  particle.  In  the  current  version  of  FastSLAM,  the  particles  are  updated  at  each 

iteration  such  that  the  carryover  of  the  previous  weights  is  not  done.  The  particle  weight  wj'”1  is 

calculated  by 


0-32 


where  the  constant  rj  which  can  be  ignored  since  the  set  of  weights  is  normalized  such  that  the 
sum  of  weights  is  unity  after  each  iteration. 

Comparison  between  EKFSLAM  and  FastSLAM  performances 

The  two  standard  Bayesian  approaches  are  explained  for  a  SLAM-based  OWLS.  It  is 
shown  that  the  EKFSLAM  linearizes  the  state  transition  and  observation  functions  and 
represents  the  joint  probability  with  a  single  high-dimensional  Gaussian;  however,  FastSLAM, 
represents  the  MN  trajectory  using  a  set  of  particles  and  conditions  the  map  on  the  MN 
trajectory.  In  particular,  the  EKF  is  subject  to  failure  where  significant  MN  uncertainty  induces 
linearization  errors.  EKFSLAM  represents  the  joint  state  posterior  PDF  by  an  approximate  mean 
and  variance.  This  can  be  problematic  from  two  aspects.  First,  due  to  linearization  error,  these 
moments  are  approximate  and  may  not  accurately  represent  the  true  first  and  second  moments. 
Secondly,  where  the  true  probability  distribution  is  non-Gaussian,  even  the  true  mean  and 
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variance  may  not  be  an  adequate  representation  of  the  PDF.  Particularly,  the  iterative  nature  of 
this  algorithm  leads  the  approximation  errors  to  accumulate  as  time  progresses  [xlii],  [xliii]. 

In  contrast  to  the  EKF,  FastSLAM  is  much  more  robust  against  linearization  error  as  the 
MN  location  can  have  an  arbitrary  PDF,  and  the  update  is  constrained  to  be  neither  Gaussian 
nor  linear.  However,  as  more  observations  are  received,  particles  which  made  poor  association 
decisions  in  the  past  tend  to  be  removed  in  the  resampling  process;  hence,  the  majority  of 
particles  tend  to  converge  onto  the  correct  set  of  associations. 

For  the  purposes  of  data  association,  FastSLAM  automatically  allows  information  to  be 
integrated  between  online  SLAM  and  full  SLAM,  as  each  particle  is  carrying  a  history  of  the  MN 
path  with  an  exponential  forgetting  memory.  Furthermore,  FastSLAM  is  simple  to  implement 
relative  to  complicated  batch  processing  algorithms.  The  important  feature  of  the  FastSLAM 
algorithm  is  that  each  particle  does  not  represent  the  current  MN  state,  but  the  entire  MN 
trajectory  and  associated  ANs  map.  This  can  have  important  effects  on  consistency  since  it 
enables  the  filter  to  accurately  estimate  uncertainty. 

However,  the  FastSLAM  particle  filter  is  really  operating  in  a  very  high-dimensional 
space:  the  space  of  the  MN  trajectories  (as  opposed  to  the  MN  location  at  a  specific  time 
instance).  Therefore,  the  required  number  of  particles  increases  exponentially  in  the  length  of 
the  trajectory.  Using  a  smaller  number  of  particles  may  result  in  underestimation  of  total 
uncertainty,  and  eventually  an  inconsistent  solution,  especially  for  the  case  whereby  the  full 
uncertainty  is  required  such  as  when  large  loops  need  to  be  closed.  Whenever  resampling  is 
performed,  for  each  annihilated  particle,  the  entire  MN  trajectory  and  map  hypothesis  is  not 
represented  in  subsequent  conditional  PDF  updates.  This  implies  that  the  random  ANs  location 
estimates  conditioned  on  these  past  locations  (particles)  are  not  sufficiently  statistically  sampled 
and  the  eventual  state  estimate  based  on  the  evolution  of  the  conditional  mean  becomes 
inaccurate.  For  remembering  long-term  uncertainty,  the  EKF-SLAM  performs  in  a  far  superior 
manner  in  comparison  with  FastSLAM  since  it  uses  a  continuous  Gaussian  representation  (the 
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mean  vector  and  covariance  matrix)  for  the  posterior  PDF,  which  does  not  degrade  as  a 
function  of  trajectory  length  [xliii],  [xlii],  However,  EKFSLAM  is  only  applicable  where  all  of  the 
PDF’s  involved  are  near  Gaussian. 

This  chapter  presented  an  overview  of  wireless  observables  from  local  wireless  sources 
which  can  be  utilized  for  an  OWLS.  Taking  advantage  of  the  signals  of  opportunity  was 
suggested  as  a  solution  for  the  indoor  environment  of  dense  urban  areas,  where  the  GPS 
cannot  provide  sufficient  accuracy.  To  this  reason,  the  deficiencies  and  problems  of 
opportunistic  wireless  localization  were  discussed.  In  summary,  the  wireless  localization 
algorithm,  which  exploits  from  an  opportunistic  signal  from  an  unregistered  local  network,  must 
combat  the  following  challenges: 

1-  The  downlink  signals  from  the  ANs  may  or  may  not  be  accurately  time  synchronized. 

2-  The  MN  clock  is  not  synchronized  with  the  ANs  clock.  Moreover,  its  clock  may  drift  in 
time  due  to  clock  instability. 

3-  Time/range  measurements  from  wireless  sources  are  subject  to  excess  delay/range 
error  due  to  multipath  fading.  Moreover,  the  range  error  is  not  stationary  due  to  the 
dynamic  feature  of  sight  conditions. 

4-  The  AN  locations  are  not  completely  known  by  the  MN. 

It  was  shown  that  the  opportunistic  wireless  localization,  in  general,  is  an  unobservable 
problem  due  to  the  inflation  of  unknown  variables,  unless  certain  conditions  can  be  assumed: 
stationary  ANs  (or  if  mobile,  they  move  in  a  deterministic  way)  and  the  MN  smooth  trajectory. 
The  first  assumption  turns  OWLS  into  an  observable  SLAM  problem  as  the  SLAM  can  take 
advantage  of  the  correlation  between  subsequent  measurements  whenever  the  same  AN  is  re- 
seen.  The  second  assumption  modeled  the  dynamic  states  including  the  MN  position,  range 
error,  and  sight  conditions  as  a  Markov  process.  The  Markovian  behavior  of  variables  can  be 
systematically  captured  by  Bayesian-based  methods.  It  was  shown  that  the  optimal  joint 
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posterior  belief  map  of  system  states  cannot  be  calculated  analytically  in  a  closed  form  due  to 
nonlinear  dynamic  nature  of  the  OWLS.  To  this  end,  two  suboptimal  Bayesian  solutions  for  the 
SLAM-based  OWLS,  EKFSLAM  and  FastSLAM,  were  introduced,  and  their  benefits  and 
deficiencies  were  discussed. 


Information  and  Estimation  Bounds  on  OWLS 

The  main  objective  of  this  chapter  is  to  derive  the  fundamental  limits  of  OWLS 
performance  in  terms  of  both  information  and  estimation  metrics.  This  chapter  illuminates 
connection  between  information  theory  and  estimation  theory  with  special  focus  on  the  SLAM- 
based  OWLS.  This  gives  an  understanding  of  what  a  specific  measurement  is  worth,  based  on 
the  metric  of  quantifying  its  information  content,  which  appears  as  an  overly  abstract  quantity 
but  is  surprisingly  practical.  For  instance,  suppose  the  issue  of  whether  a  MEMS  IMU  can  be 
added  to  the  MN  receiver.  This  will  provide  additional  information  regarding  the  update  of  the 
position  of  the  MN  in  the  SLAM  algorithm.  Is  this  worthwhile?  How  much  additional  information 
does  the  IMU  provide?  A  component  of  this  thesis  is  the  consideration  of  information  in  regards 
to  the  constraints  of  the  SLAM  formulation  in  the  context  of  the  OWLS. 

In  estimation  theory,  the  MMSE  and  CRLB  are  considered  as  the  most  fundamental 
estimation  measures  to  illustrate  how  accurately  each  individual  parameter  can  be  recovered 
from  channel  observables.  Estimation  theory  is  tightly  coupled  with  the  averaged  squared  error 
of  estimation,  regardless  of  whether  the  unknown  state  q  is  a  deterministic  parameter  (as  in 
classical  approaches),  or  whether  the  unknown  state  q  is  a  random  variable  (as  in  Bayesian 
approaches)  [xliv],  [xlv],  given  by 

Unbiased  Classical  Estimator  :  var(q(z))  =  £[[q(z)  -  qj2  J  0-33 

Bayesian  Estimator  :  MSE(q(z))  =  2?  [[ q(z)-q]2] 
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However,  to  unify  the  concept  of  estimation  error  in  classical  and  Bayesian  philosophies, 
we  use  MSE  as  a  general  term  for  the  averaged  squared  error  of  estimation  for  both  classical 
and  Bayesian  methods,  accepting  the  difference  that  the  averaging  is  with  respect  to 
observation  PDF  (since  there  is  no  prior  information  available  about  non-random  parameter)  in 
classical  philosophy,  while  the  averaging  is  with  respect  to  joint  PDF  of  observations  and 
random  parameter  in  Bayesian  philosophy. 

Cramer  Rao  Lower  Bound  (CRLB)  theorem  says  that  the  MSE  corresponding  to  a 

unbiased  estimator,  q(z) ,  of  a  non-random  parameter  q,  given  the  observation  vector  z  ,  can 
be  lower  bounded  [xlv],  [xlvi],  as  follows: 

£{(q(z)-q)(q(z)-q)r}>J-1 

where  J  is  the  FIM  of  non-random  parameter  q  ,  obtained  as 
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J  D  [  Vq  ln  P z|q  (Z  I  q)]  [  Vq  ln  P z|q  (Z  I  q)] 


However,  the  extended  CRLB  theorem  for  random  parameters  (also  known  as  Bayesian 
CRLB,  posterior  CRLB,  and  Van  Trees  inequality)  indicates  that  the  MSE  of  an  unbiased 
estimator  q(z)  can  also  be  lower-bounded  [xlvi],  [xlvii],  as  follows, 

0-36 


£{(q-q)(q-q)  } ^ -Kl 

where  J tot  is  total  Bayesian  FIM  (BFIM)  of  the  random  parameter  q  ,  defined  by 


J  tot  ^  J Z  +  J P 
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and  Jz  is  the  measurement  information  matrix,  representing  the  information  data  obtained  from 
the  observables,  z, ,  stated  as 

in  0-38 


JZ  D  Kq  [  Vq  ln  Pz\q  (Z  I  q)][Vq  ln  P,\q  (z  I  q)] 


The  expectation  is  taken  under  the  joint  PDF  of  measurements  and  states.  The  matrix 
J p  denotes  the  apriori  information  about  q, ,  whose  elements  are  defined  as 
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4  34 


d^Pq (q)  d In pq (q) 

-  p 

d2  In  A?(q) 

dqt  Gqf 

dq,q, 
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,  i,j  =  \,..K 


where  K  =  dim(q,)  stands  for  the  state  space  dimension.  It  must  be  noted  that  Bayesian 

CRLB  theorem  holds  under  some  mild  regularity  conditions,  implying  that  the  joint  PDF  does  not 
have  infinite  moments. 

Information  theory  provides  a  deep  analytical  approach  to  assess  tracking  performance 
based  on  the  available  information.  Information  theory  application  in  tracking  problems  is 
traditionally  known  by  Fisher  information  (FI)  regarding  the  CRLB  calculation.  Fisher  information 
is  recognized  as  a  way  to  measure  the  amount  of  information  that  an  observable  carries  about 
the  parameters  to  be  estimated  upon  either  their  joint  (for  random  parameters)  or  likelihood  (for 
non-random  parameters)  probability  density  functions. 

Another  important  metric  in  information  theory  is  the  input-output  mutual  information 
(Ml)  that  is  used  as  an  indicator  of  how  much  coded  information  can  be  pumped  through  a 
channel  reliably  given  a  certain  input  signaling.  However,  the  main  application  of  Ml  has  been 
generally  known  in  coding  operations,  Guo’s  theorem  [xi]  shows  that  Ml  can  also  represents  a 
powerful  measure  to  provide  a  deep  insight  into  the  performance  of  estimation  problems  such 
as  the  SLAM.  Mutual  information  I(X,Y)  can  be  considered  as  a  global  measure  of  available 
information  in  channel  observable  Y,  about  state  variable  to  be  estimated  X  ,  since  it  does  not 
make  any  assumptions  regarding  the  receiver  processing  of  Y .  However,  the  minimum  MSE 
(MMSE)  criterion  is  obtained  according  to  the  quadratic  cost  function  of  error;  this  implies  that 
posterior  probability  function  needs  to  be  of  an  exponential  (Gaussian)  type  to  reveal  all  related 
information  about  estimation  error  from  the  second  order  moment.  As  a  consequence,  CRLB  or 
BCRLB  may  not  fully  characterize  the  accuracy  of  localization  from  a  non-Gaussian  posterior 
PDF  which  requires  higher  order  moments,  in  addition  to  the  mean  and  the  covariance. 
Nevertheless,  BCRLB/CRLB  and  BFIM/FIM  is  regarded  as  a  useful  metric  not  only  to  assess 
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the  performance  of  implemented  algorithms  but  also  to  predict  the  best  achievable  performance 
as  a  system  design  tool. 

This  chapter  introduces  both  the  information  and  the  estimation  metrics  and  their 
relations  to  evaluate  OWLS  performance  in  terms  of  efficiency,  consistency  and  observability. 
Information  matrix  is  introduced  as  a  quantitative  and  visual  benchmark  for  all  the  disparate 
sources  of  information  utilized  in  OWLS.  In  this  regard,  the  contribution  of  every  source  of 
information  in  a  SLAM-based  OWLS  is  explained  by  some  basic  cases.  Based  on  the  system 
parameters  introduced  for  a  range-based  OWLS  in  Chapter  2,  four  system  models  are  proposed 
to  present  a  complete  set  of  system  state.  BFIM  is  derived  for  each  proposed  system  model  as 
a  fundamental  limit  of  accuracy  for  a  range-based  OWLS.  Later,  the  relation  between  FI  and 
other  estimation  metrics  is  explained  which  can  be  exploited  as  useful  tools  to  assess  the  SLAM 
algorithms  in  terms  of  consistency  and  localization  accuracy.  Finally,  an  applicable  example  of 
Ml  for  a  SLAM-based  OWLS  is  provided  through  introduction  of  insightful  connection  between 
Ml  in  information  theory  and  MMSE  in  estimation  theory  from  Guo’s  theorem. 


Fisher  information  and  BCRLB  application  in  a  SLAM  problem 

The  concept  of  the  Bayesian  information  matrix,  denoted  by  J  ,  is  key  to  SLAM 
algorithms.  The  SLAM  filter  efficiency  can  be  investigated  by  comparing  the  extracted 
information  with  the  existing  information,  quantified  as  J .  It  is  loosely  regarded  as  a 
quantitative  summary  of  all  information  regarding  the  unknowns  in  one  matrix.  The  relation  to 
quantifiable  information  is  only  valid  for  the  case  where  the  problem  is  linear  and  jointly 
Gaussian  (LJG);  however,  for  problems  where  the  nonlinearity  is  mild  and  the  measurement 
and  state  update  noise  are  approximately  Gaussian,  it  is  still  applicable. 

The  MLE  of  the  unknown  variable  vector  q  can  be  typically  approximated  as  being 

jointly  Gaussian  represented  by  w(q,  J .  Here,  q  represents  the  actual  value  of  the  state 
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variable  vector.  For  typical  problems  that  are  not  jointly  Gaussian,  the  PDF  approximation  of  the 
MLE  as  7V(q,  J _1)  is  still  generally  valid.  Regardless,  J  is  the  information  matrix  and  is  the 
inverse  of  the  covariance  matrix  of  the  MLE  estimators  for  q .  In  the  OWLS,  whenever  the  LJG 
assumption  is  justified,  the  Fisher  information  matrix  can  be  derived  as  the  fundamental  limit  of 
localization  and  mapping  accuracy  for  the  OWLS  problem. 

Information  sources  for  a  SLAM  problem  is  generally  derived  from  the  MN  motion 
process,  measurements  from  ANs,  and  the  prior  knowledge  about  the  MN  initial  state  and  AN 
locations.  To  understand  the  contribution  of  each  information  source  on  the  SLAM  problem,  the 
FIM  calculation  process  is  explained  by  6  simple  linear  Gaussian  cases  while  the  problem 
complexity  is  eventually  increased  for  next  case. 

In  the  first  case,  consider  the  MN  undergoes  random  motion  in  the  x  direction  according 
to  the  first-order  Markov  (random  walk)  model.  Assuming  that  the  initial  state  is  distributed  as 
x0  □  7V(o,(j02),  that  represents  the  prior  PDF.  The  update  step  from  x0  to  xl  is  illustrated  in 
Figure  0-13 


Figure  0-13  One-step  state  update  in  a  ID  linear  problem 


The  update  model  is  defined  by 


0-40 

x1  =  x0  +  u 

where  the  process  noise  is  distributed  by  u  □  ./v(o,cr2) .  The  joint  posterior  PDF  of  {x0,xj  is 

given  by  a  normal  PDF  with  the  zero  mean  and  the  covariance  matrix  of 

r<72  C72  1  °-41 

u0  u0 

^  ~  ^2  2  ,  2 

[(Jo  CT0  +  Gu  J 

Note  that  the  covariance  of  x1  given  x0is  denoted  as  cr2  ,  determined  by 
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_2  _/  2, ,_2\  ^0  _  2 

CTI|0  K  +  <T,  )  2 


Likewise,  the  covariance  of  x0 1  jCj  is  given  by 


_2  _  -.2 
^011  <T> 


CTn 


_2  .  _2 
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2  2 

^0  +  <*„ 
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which  for  the  case  where  cr2  □  cr2  ,  is  approximately  cr2 .  If  there  is  no  measurement  then  the 
variance  of  the  estimate  of  x1  is  cr2  +  cr2 .  The  Fisher  information  matrix  is  obtained  by 


J  =  QT'  = 


2  ,  _2 
°~0  +Cr„ 

2  2 

Wu 


-2 


-cr 


-2 
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—cr  <j 

u  u 

It  is  notable  that  diagonal  elements  of  the  information  matrix  are  equal  with  the  inverse  of 


conditional  covariances  as  follows, 


[J]u=< 

[42=<v 

Now,  in  the  second  case,  consider  another  problem  where  only  the  measurement 
associated  with  the  anchor  point,  ml ,  is  available,  as  illustrated  in 

Figure  0-14. 


0-45 


Figure  0-14  Single  measurement  at  a  single  point  trajectory 


With  the  assumption  of  z,  □  7V(x0,cr),  the  covariance  of  x0  given  the  measurement  is  stated 
as 
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In  the  third  case,  consider  the  measurement  as  an  amendment  to  the  prior  one  step 
update  problem,  as  shown  in  Figure  0-15.  Similarly,  the  covariance  of  x0 ,  given  the  previous 


state  and  measurement,  is  computed  by 


<T, 


1|0  ,z 


=  crn,  +<r„  = 


2 

a;+a- 


+  cr~ 


0-47 


Figure  0-15  Single  update  step  with  a  single  measurement 


The  FI  calculation  procedure  can  be  viewed  differently  with  the  procedure  based  on  the 

state  vector  of  q  =  [x0  x,]r  .  Calculation  of  total  information  involved  for  the  above  cases  can 

be  determined  by  the  following  steps: 

1 .  First,  start  with  the  covariance  matrix  of  the  state  update: 


Q  = 


2  2  2 

_cr0  cr0  +  cr, 
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2.  The  information  matrix  from  the  measurement  is  given  by 

°  ° 

as  it  refers  to  the  x0  position  state. 

3.  Compute  the  total  information  as  the  sum  of  prior  information  Jp  =  Q  1  and 
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measurement  information  J_,  expressed  by 
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J tot  J P  + 


a;2  0 

0  °_ 

4.  Determine  the  inverse 


Qtot  ~  J tot 

For  the  current  case,  the  posterior  PDF  covariance  matrix  is  determined  as 
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2 _ 2 


Qtot  = 


Uot  ~ - 1  ^-2 

cr0  +crz 
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Note  that  the  second  diagonal  term  of  Qtot  defines  the  estimation  variance  of  x1  given  x0  and 


z  as  : 


[  Qt»t  ]22 


2.  2.  ,  2  2  ,  2  2 

cr0crr  +o-zo-„  +a-0o-„ 

_2  ,  2 
+(Tz 


_.2  2 
o-02  +  a. 


,  2  2 
2  +  ~  °"l|0,z 
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Incrementing  the  information  matrix,  J  ,  with  a  new  incremental  entity  of  information 


such  as  a  measurement,  is  shown  in  this  example.  The  only  processing  problem  is  that  J  has 


to  be  eventually  inverted  to  extract  the  actual  variances  of  the  state  variables.  Of  course,  the 
other  issue  is  that  of  the  underlying  assumption  of  the  problem  being  jointly  Gaussian.  Note  that 


if  the  measurement  provides  no  information,  i.e.  cr2  — >  oo  ,  then  crj|0  _  — >  cr(,  +cr2  as  in  the 


previous  example. 

Now,  in  fourth  case,  it  is  assumed  that  the  AP  independent  measurements  are  available 
from  both  x0  and  x, ,  as  illustrated  in  Figure  0-16. 


Figure  0-16  Two  independent  measurements  of  two  subsequent  states 
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The  information  matrix  for  the  two  independent  measurements  is  given  as 


J  = 


-2 


0-54 


The  total  covariance  can  be  computed  similarly  as 


Qtot  ~  J tot  ~ 


( 
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In  an  estimation  problem,  a  parameter  is  said  to  be  observable  if  there  is  sufficient 
information  to  estimate  it  without  ambiguity.  FIM  can  also  be  used  to  investigate  system 
observability.  In  order  to  have  an  observable  system,  the  FIM  must  be  invertible.  If  the  FIM  (or 
BFIM)  is  singular  or  very  ill-conditioned  then  CRLB  (or  BCRLB)  does  not  meaningfully  exist. 
As  shown  in  Eq.(0-48),  the  single  observation  of  an  AN  results  in  a  singular  FI  matrix;  it  means 
that  the  amount  of  information  is  insufficient  for  the  estimation  problem.  However, 
unobservability  can  be  resolved  if  the  MN  moves,  as  shown  in  Eq.  0-49,  and/or  the  AN  is 
observed  in  more  than  a  single  MN  state,  as  shown  in  Eq.  (0-53). 

In  the  fifth  case,  the  previous  localization  problem  is  transformed  to  the  SLAM  problem 
where  the  measurement  is  obtained  from  a  feature  point  with  initially  unknown  location,  as 

illustrated  in  Figure  0-16.  The  state  vector  is  modified  as  q  =  [x0  x,  m  f ,  where  the 

covariance  matrix  of  the  updates  is  given  by 
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The  information  added  by  the  measurement  between  m  and  x0  is  captured  in  the 


information  matrix  as 


X  o  A 


0  0 


0 


~  °  -4 

a:  a: 
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where  the  component  of  the  matrix  is  calculated  by 


\J^=~E 


In  (jV(/w-x0,<7_:)) 
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In  a  /similar  way,  the  information  added  by  the  measurement  between  m  and  xl  is  given 


by 


J,.  i  = 


0  0 
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Hence,  the  overall  covariance  matrix  is  obtained  by 


Qtnt  = 
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An  interesting  point  in  this  example  is  that  however,  the  observation  information  matrix  of 
each  state  Jzi,i  =  1, 2 ,  is  a  singular  matrix,  the  multiple  observation  of  a  FP  in  a  correlated 

trajectory  makes  the  estimation  problem  observable.  It  must  be  noted  that  a  single  FP 
measurement  that  is  only  connected  to  one  state  variable  does  not  convey  information. 

By  now,  a  systematic  approach  to  calculate  information  connectivity  of  FPs  and  the 
motion  update  from  one  time  step  to  the  next  is  established.  It  is  shown  that  the  information 
conveyed  is  therefore  minimal  if  each  FP  is  assumed  to  have  only  two  associated  observations; 
hence,  the  information  matrix  is  sparse. 

The  final  case  complicates  the  SLAM  problem  to  bring  up  the  localization  problem  in 
cooperative  networks  where  the  mobile  nodes  can  help  each  other  by  sharing  their  information. 
Although  cooperative  SLAM  problems  are  out  of  scope  of  the  current  research,  a  simple 
illustrative  example  is  given,  as  illustrated  in  Figure  0-17.  The  FIM  calculation  procedure,  as 
systematic  way  of  accounting  for  the  diverse  set  of  constraints,  is  extended  for  two  MN 
cooperative  tracking  problem  which  two  MNs  meet  a  single  FP. 


Figure  0-17  Two  cooperative  MN  localization  problem 
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The  update  covariance  matrix  is  calculated  by 
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Similar  to  previous  example,  the  measurement  information  matrix  is  given  by 
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The  total  covariance  matrix  and  so  the  information  is  computed  as  follows 
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The  diagonal  elements  of  Qtot  are  defined  by 
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In  a  dense  obstacle  environment,  the  opportunistic  sources  may  not  be  able  to  provide 
sufficient  localization  information  to  the  MN  because  of  radio  blockage  and  limited  range.  In 
such  situations,  cooperation  between  MNs  can  significantly  improve  opportunistic  localization 
performance.  Cooperative  OWLS  can  be  investigated  by  the  case  where  the  FPs  are  also 
moving  and  have  their  own  update  motion  models,  as  shown  in  Figure  0-18.  In  this  case,  the 
new  states  associated  with  FPs  will  be  augmented  to  the  state  space  as 


y=[> 


x, 


mn 
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X  .  For  each  time  step,  the  covariance  matrix  for  the  update  state  will  be 
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Figure  0-18  Cooperative  SLAM  with  moving  FPs 
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In  similar  way,  the  measurement  information  matrix  is  shown  as 
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A  systematic  approach  to  calculate  FI  elements  and  unify  all  information  from  different 
sources  of  information  into  one  single  matrix  is  explained  with  some  simple  linear  Gaussian 
cases.  In  the  next  subsection,  the  FI  concept  is  generalized  for  the  nonlinear  OWLS  problem  in 
a  SLAM  framework. 


Opportunistic  wireless  localization  system  models 

Consider  a  single  MN  moving  smoothly  through  a  2D  space  in  a  random  trajectory  and 
receiving  signals  from  an  opportunistic  wireless  network  that  consists  of  an  arbitrary  number, 
Nan  ,  of  stationary  ANs.  The  location  of  ANs  may  not  be  known  to  the  MN;  however,  their 

correspondence  is  not  an  issue;  therefore,  the  additional  mapping  variables  associated  with 
mapping  do  not  have  to  be  carried  in  the  state  vector.  Figure  0-19  shows  a  2D  illustration  of  an 
OWLS  scenario,  where  blue  lines  represent  the  continuous  MN  trajectory,  and  red  circles  on  the 
MN  trajectory  represent  the  MN  locations  where  the  measurements  are  taken  and  state 
estimates  are  updated.  FPs  and  APs  are  shown  by  red  and  green  squares,  respectively,  and 
LOS  blockage  is  created  by  a  wall  shown  by  a  light  green  rectangular. 

The  MN  location  p,  is  modeled  as  a  random  walk  process  ruled  by  the  system  equation 
p,  =  p(_,  +  v,,  where  the  v,  e  □  2  is  the  motion  driving  process  with  the  known  zero  mean 
Gaussian  distribution,  fv{yt) ,  or  v,  □  N(0,Qp ) .  The  first-order  Markov  model  for  motion  is  the 
most  non-informative  model  for  the  smooth  MN  trajectory  assumption,  which  considers  the 
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velocity  to  be  zero  and  fv(yt)  to  describe  the  velocity  variation.  Higher-order  motion 

parameters  such  as  velocity  and  acceleration  can  also  be  modeled  for  the  MN  motion  process 
for  a  smoother  and  more  informative  trajectory  assumption.  Moreover,  the  motion  control  data 
or  trajectory  estimation  from  CV  or  IMU  sources  can  be  replaced  by  the  random  walk  model  to 
insert  more  information  in  tracking  algorithms,  whenever  available. 


Figure  0-19  MN  tracking  in  mixed  LOS/NLOS  conditions 

Let  n4P  denote  the  number  of  all  APs  whose  locations,  represented  by 
{ mjx , mi  y  | ,  /'  =  1 :  NAP ,  are  precisely  known  to  the  MN  ,  and  nfp  denote  the  number  of  FPs 

whose  locations  {mix,miyj,i  =  NAP  +  1 : NFP  +NAP  are  unknown  to  the  MN.  The  range 

measurements  obtained  from  opportunistic  signal  reception  from  ANs  can  be  subject  to  range 
offset  error  due  to  non-synchronization,  multipath,  and  NLOS/LOS  effects. 

In  this  research,  the  SLAM-based  opportunistic  localization  problem  is  investigated  in 
four  system  models.  In  the  first  case,  it  is  assumed  that  the  MN  receives  a  synchronized  signal 
from  a  synchronized  network  of  ANs.  In  addition,  ANs  provide  LOS  measurements,  or  if  not, 
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ANs  sight  conditions  are  known  to  the  MN,  and  NLOS  measurements  are  discarded  from  the 
measurement  data  set. 

The  pseudo  range  measurement,  obtained  at  the  ?th  time  instant  from  the  i  th  AN,  is 
explicitly  presented  by 

z.  =  Ip  ,-m.l  +  w. ,  0-67 

i,t  t  i  |  l,t 


where  p,=[x,,v,]/  is  the  MN  location  vector  at  time  step  t,  and  m  i=\_mix,mijj  is  the  z'-th 
AN  location  vector.  wit  is  an  additive  white  Gaussian  noise  with  a  mean  power  i?  jw,.,2  j  =  cr2,. . 
In  this  case,  the  state  vector  is  denoted  by  q,  =  [p,  J  .  where 


mfP  = 


mA 


m 


Njp+i  ,y 


mA 


m 


NAN,y 


is  the  stacked  vector  of  the  FP  locations.  Since 


the  ANs  are  stationary,  the  transition  PDF  from  the  MN  position-  ANs  position  state 

q?_t  =  [p,_i  m//;]  to  the  next  one  q,  =  [p,  m/p]  is  defined  only  by  the  MN  motion  update, 


p(&Mt-x)  =  p(v,  Ipm)- 

The  second  system  model  assumes  that  synchronized  ANs  are  in  the  MN’s  line  of  sight 
but  are  not  synchronized  with  the  MN.  The  range  measurement  obtained  at  the  Mh  time  step, 
from  the  i  -th  AN,  is  explicitly  presented  by  the  pseudo-range  biased  due  to  the  random  MN-AN 
clock  offset  as 

zu  =  ~ -T )2  +  {mi.y  ~ T, )2  +  bt  +  Wut  0  68 

i  1, 2,  ...5  N ^ 

where  bt  is  the  random  range  offset  which  takes  into  account  both  the  deterministic  clock  time 

offset  between  the  ANs  and  MN  and  random  MN  receiver  clock  instability  properties.  The 
random  range  offset  is  modeled  by  an  independent  random  walk  from  the  MN  motion  process 
as  p(bt  1 N(0,al)  where  cr2  is  commensurate  with  Allen  variance  of  the  MN  clock.  In  this 
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case,  the  state  vector  is  defined  by  q,  =  [p,  bt  m/p  J  .  If  the  ANs  are  not  synchronized,  the 
bias  variable  for  each  AN  will  be  different  and  the  state  vector  is  modified  by 


-\T 


p,  b,  mfp  ,  where  b,= 


u 


NanJ 


be  studied  as  the  third  system  model,  given  by 


is  the  vector  of  ANs  range  offsets;  this  will 


0-69 


/( mi.x  -  x,  V  +  ( mi.v  -  y,  V  +  bu  +  wu 

i  =  1,  2,...,NaJ~^ 

The  clock  drift  process  of  each  AN  evolves  independently  from  the  clock  drift  processes 


of  other  ANs,  so  the  transition  PDF  of  range  offset  variables  is  given  by  p( b,  |b  AT0,&), 


where  Qh  =  diag  |  ai  cr, 


The  fourth  system  model  considers  a  general  multipath  scenario  where  the 
measurement  is  also  affected  by  non-stationary  channel  behavior.  Multipath  and  LOS  blockage 
are  the  most  penalizing  error  sources  in  wireless  localization.  It  can  be  classically  modeled  by  a 
white  noise  in  the  signal  bandwidth.  NLOS  measurements  are  detected  as  an  abrupt  change  in 
mean  and/or  variance  of  a  Gaussian  noise  [xlviii].  However,  if  the  NLOS  condition  is  left 
undetected,  it  substantially  degrades  localization  accuracy.  One  approach  to  deal  with  the 
NLOS  conditions  is  to  introduce  a  noise  model  for  LOS  propagation  and  a  noise  model  for 
NLOS  propagation,  where  the  transition  between  the  LOS  and  NLOS  mode  is  modeled  with  a 
Markov  chain.  In  this  case,  the  state  vector  is  characterized  by  the  MN  position,  FP  locations, 


range  offset  vector,  and  ANs’  sight  conditions  s. 


,  given  by 


q,=[p,  b,  s,  m  f 

where  each  AN  sight  condition,  t ,  is  modeled  by  a  binary  Markov  Chain,  with  the  value  sit  =  0 


is  assigned  to  the  LOS  event,  and  sit  =1  is  assigned  to  the  NLOS  event.  As  earlier  discussed, 
in  mixed  NLOS/LOS  conditions,  the  suggested  range  measurement  model  is  characterized  by 
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introducing  sight  state-dependent  measurement  noise  distributions  for  fourth  system  model, 
given  by 


zu,  =  J ( mL*  -x,)2  +(mLr-yt)2  +  bu  +  wu 
i  =  1,2 

where 

n  \mK J  <fsu  =0  0-72 

\N^i,nlos^i,nloS)  Si,t=l 

Under  LOS  condition,  the  measurement  is  only  corrupted  by  additive  receiver  noise, 
which  is  described  by  a  zero-mean  Gaussian  distribution  with  a  variance  of  aflos .  The  NLOS 
excess  delay  is  modeled  by  the  abrupt  jump  in  the  mean  and  variance  of  measurement  noise, 
where  o]nlm  >  crfJos  and  fii  nlos  >  0  [xlix],  [I],  [li].  We  assume  that  the  M„nios ,  °lnlos  ,  and  afhs  are 

known  from  prior  information  about  the  room  size  and  the  calibration  stage. 

The  Nan  Markov  chains  are  combined  into  a  single,  augmented  Markov  chain,  denoted 

by  s,  which  have  2Nan  possible  states.  As  long  as  the  range  measurements  are  collected  from 
ANs  at  different  locations,  the  LOS/NLOS  transition  can  be  assumed  to  be  independent. 
Therefore,  the  transition  probability  matrix  FI  of  the  augmented  Markov  chain  can  be  expressed 
in  terms  of  the  transition  probability  matrices  of  n,-,/  =1,...,#^  of  each  AN,  according  to 

n  =  ni®n2®...®n  v,v .  Where  0  denotes  the  Kronecker  product.  The  description  of  four 
proposed  system  models  is  summarized  in  Table  0-1. 


Table  0-1  System  models’  description 


System  Model 
scenario 

ANs  location 

MN  clock  instability 

ANs 

synchronization 

NLOS  conditions 

i 

Unknown 

The  MN  are  synchronized 
with  ANs 

ANs  are 
synchronized 

Sight  conditions  are  known 
and  NLOS  measurements 
are  discarded 

Solution:  ANs  unknown 
locations  are  treated  as 
FP’s  of  unknown 
location. 

2 

Unknown 

The  MN  are  not  synchronized 
with  ANs 

ANs  are 
synchronized 

Sight  conditions  are  known 
and  NLOS  measurements 
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Solution:  similar  to  the 
system  model  1 

Solution:  MN  clock 
instability  is  treated  as 
random  first-order  Markov 
process. 

are  discarded 

3 

Unknown 

The  MN  are  not  synchronized 
with  ANs 

ANs  are  not 
synchronized 

Sight  conditions  are  known 
and  NLOS  measurements 
are  discarded 

Solution:  similar  to  the 
system  model  1 

Solution  For  every  MN-AN  link,  the  clock  drift  is 
treated  as  a  random  first-order  Markov  process 

4 

Unknown 

The  MN  are  not  synchronized 
with  ANs 

ANs  are  not 
synchronized 

Sight  conditions  are 
unknown 

Solution:  similar  to  the 
system  model  1 

Solution:  similar  to  the  system  model  3 

Solution:  Sight  conditions 
are  modeled  by  Boolean 
switch. 

Fisher  information  for  a  SLAM-based  OWLS 


A  direct  approach  to  quantify  the  potential  accuracy  of  a  trajectory  estimation  based  on 
the  observed  data  samples  is  formulated  by  the  Fisher  information  matrix.  However,  as  stated 
earlier,  the  quantifiable  information  is  only  valid  where  the  problem  is  linear  and  jointly  Gaussian 
or  at  least  nonlinearity  is  mild  and  the  measurement  and  state  update  noise  are  approximately 
Gaussian.  FI  determines  a  fundamental  limit  for  localization  accuracy  by  using  information 
inequality  through  the  CRLB  theorem.  This  information  bound  predicts  the  best  achievable 
performance  even  before  the  OWLS  is  built,  so  that  it  can  be  utilized  as  a  system  design  tool. 

For  the  SLAM-based  OWLS,  where  q0.,  represents  the  history  of  unknown  system  state 

vector,  the  BCRLB  as  a  lower  bound  on  the  MSE  matrix  of  any  unbiased  estimator  q0.,(zl7)  is 
obtained  by  substituting  q0.(  by  q  and  zw  by  z  in  Eq.  (0-35)-(0-38),  yielding 

E  {(q0:,  -  q0:, )  (q0:,  -  )r }  ^  -4'  (q(w ) 

where  Jtot( q07)  denotes  the  BFIM  in  regards  with  history  of  state  vector.  Jtot{ q07)  is  quantifiable 
information  associated  with  full  SLAM  processing;  however,  in  OWLS,  the  current  system  state 
q,  is  of  particular  interest.  It  can  be  easily  proved  [lii]  that  BCRLB  for  q,  is  given  by  the  KxK 
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lower  right  submatrix  of  Jt0]( q0t),  where  K  =  dim(q,) .  Let  decompose  q,  as  , 

«/to,(q0:J  can  be  written  by 


J tot  (q0:r ) 


A  b, 
Bj  C,_ 


0-74 


\  E{-^^p*,Az i*.qte)} 

E{-A^-'\npzq(zl.t,q0.,) 

[  ^{-Aq;lnPZ,q(Zl,>qO:,)} 

where  Af  is  the  operator  of  the  second  order  derivative  and  is  also  defined  versus  the  first- 
order  gradient  as  Af  =  .  Given  Eq.  (0-73)  decomposition,  the  Fisher  information  matrix 

regarding  the  current  state  q, ,  Jtot( q,) ,  is  obtained  by 


AoMl)  =  Ct-BtTAt-lBt 

Generally,  there  are  two  types  of  algorithms  that  calculate  the  Fisher  information  for 
online  processing.  The  first  type  avoids  inversion  of  the  large  matrix  of  Jlot{ q0v)  and  calculates 


the  Jtol( q,)  recursively  based  on  Eq.  (0-93)  [xii],  [lii]  as 

0-76 

J,Mk)  □  Df  -  Dl\D"  +  4(qH)r‘^ 

where 

^(q*)D  D?-D?(D\'+JtoX^)YlD?  0-77 

L>?  =E{-Al\np(qk+l\qk)} 

Dl2  =  E  {-Al+' In  p( q*+Jq*)} 

A21=^{-A^lnMq,+1|q,)}  =  [^2f 

A22  =£{-Aq!i;  ln^(q*+i  I  q*)}+£{-Aq“  lnP(zk+ 1 1  q*+i)} 

For  a  LJG  case,  the  recursion  for  Jtol( q,)  is  identical  to  inverse  of  the  posterior  PDF 

covariance  matrix  calculated  from  the  Kalman  filter.  However,  in  a  nonlinear  SLAM  problem, 
due  to  correlation  between  state  vector  estimates,  as  earlier  discussed,  the  recursive  methods 
partially  discard  the  correlated  information  between  the  current  state  and  older  states. 
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The  second  types  of  methods  focus  on  calculating  full  SLAM  BFIM  Jtot(%.t)  numerically 


[Iv]  or  based  on  reasonable  approximations  [liv],  For  the  problem  at  hand,  we  utilized  the  second 
approach  since  the  inversion  of  Jtot{ q07)  is  not  an  issue  for  the  dimension  of  the  simulated 

problem;  otherwise  the  recursive  methods  need  to  be  used.  Except  the  forth  system  model 
which  is  described  by  a  hybrid  discrete  and  continuous  state  space  and  must  be  treated 
differently,  the  proposed  system  models  for  the  SLAM-based  OWLS  can  be  described  by  a 
linear  model  for  state  update  and  a  nonlinear  model  for  the  measurement  function,  yielding 

q,=q  °-7S 

zu,=hM  +  wi,>  »  i  =  h2,...,NAN 

where 


0-79 


F 


1  0  0  0 

0  1  "  • .  :  0 

:  •.  •.  o  : 

0  ...  0  1  0 


0 

0 

0 


Ld  and  Ls  denote  the  dimension  of  dynamic  variables  and  stationary  variables  of  state 
vector,  respectively.  The  noise  vectors  are  mutually  independent  white  process  distributed  as 
\d  ,□  N(0,Qd)  and  wit  □  7V(0,cr2_).  Qd  is  a  nonsingular  covariance  matrix  of  joint  dynamic 


variables  process  where  for  the  system  model  1  ,  Qd  -  Qp ;  for  the  system  model  2, 

Qd  =  diag {0p,cr62}  and  for  the  system  model  3  and  4,  Qd  =  diag\Qp,Qb  j .  It  is  assumed  that 
the  MN  clock  instability  process  evolves  independently  from  the  MN  motion  process.  From  Eq. 
(0-91),  the  BFIM  </to((q0.,)  for  a  SLAM-based  OWLS  is  given  by 

,  ,  0-80 

J,0  (q<*)  =  E  {-Aq“  in  a,  ,(zto.qw)) 

According  to  Eq.  (0-35)-(0-38),  the  total  BFIM  is  decomposed  using  Bayes’  rule,  stated 
as 
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J tot  (*lo:r)  J/  07  )  ^  J ) 


0-81 


where  Jz( q0.,)  denotes  the  BFIM  of  measurement  data,  defined  as 


0-82 


and  where  /P(q0.,)  denotes  the  BFIM  of  prior  data,  defined  as 


*4(qto)  =  {-A^:  In  pq(q0:,)} 


0-83 


The  evaluation  of  BFIM  of  the  measurement  data  involved  the  computation  of 
expectation,  which  is  difficult  to  express  analytically  for  the  nonlinear  model  given  in  Eq.  (0-77) 
unless  some  reasonable  approximations  can  be  assumed.  We  utilize  first-order  Taylor 
expansion  to  construct  a  linear  approximation  of  h:(.)  around  q,  which  is  the  most  likely  system 
state  deemed  at  the  time  of  linearization: 


0-84 


3q, 


liberalized  measurement  process  is  given  by 


0-85 


Using  the  linearized  measurement  model  approximation,  the  BFIM  from  measurement 
data  can  be  expressed  by 
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0-86 


Jz(<\0,)  =  iE {-A^“  In p( zk  I  q, )} 

k= 1 

k= 1 

k= 1 

where  <^’7,  denotes  a  (£  + 1) x  (k  + 1)  dimensional  matrix  whose  elements  are  all  zero  except  at 

the  i -th  row  and  the  y'-th  column  which  is  one,  and  Jz( qk)  denotes  the  BFIM  portion  from  the 

measurement  data  at  the  A:-th  time  step.  It  must  be  noted  that  the  FPs  are  stationary,  so  unlike 
dynamic  state  variables  such  as  MN  location  variables,  there  is  no  need  to  define  a  new  state 
vector  portion  regarding  to  FP  location  variables  at  each  time  step.  To  unify  the  information 
regarding  FP  location  variables  from  all  time  steps,  the  state  vector  variables  are  divided  in  two 

dynamic  and  stationary  portions  as  =  [qrf4 :,qs.  J^and  the  history  of  state  vectors  is  re-sorted 

as  q0/t  =  [qt/(H.,qv]r  ■  In  this  regard,  the  BFIM  for  mixed  stationary  and  dynamic  states  from 
measurement  data  is  modified  by 

•4 too*)  =  Z^{-Aq“  ln^(z *1  <!*)} 

k= 1 

=  X£{-C2U*‘®^"lnP(Zl|q,)} 

k= 1 

+  £{-<5,™*'®^"  lnp(zt  |  q,)} 

+  E{-S;T2  ®a;;  In  p(z,  |  q„)} 

To  calculate  the  prior  information,  we  need  to  consider  the  update  model  of  dynamic 
state  variables  as  well  as  initial  conditions.  Since  the  FPs  are  assumed  to  be  stationary,  the 
update  information  is  only  applied  for  dynamic  variables.  The  initial  conditions  are  divided  into 
two  portions:  the  stationary  portion  which  is  later  added  to  the  measurement  information  to 
calculate  the  total  information,  and  the  dynamic  portion  which  is  combined  with  the  update 
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information.  According  to  the  LJG  model  of  the  state  update  process,  and  assuming  that  initial 


state  q0  is  Gaussian  distributed  as  q0 1,.Y(0.0>,)  which  is  composed  of  initial  condition  for 
dynamic  variables,  q0^  0  N(0,Qd0 ) ,  and  stationary  variables,  q0i  □  N(0,Qs0),  the  prior 
information  about  dynamic  states  q</0  and  qrfl  with  the  update  state  model 

q</,i-q</,oD  N(Q,Qd)  is  given  by 

•M<uM>=a(q,«,r‘  °'8! 


Qu(  q 

d,  0:1 ) 

7 

d 

l _ 

o" 

+Jd 

0 

0 

where 

Jd  = 

d 

1 _ 

i 

P 

i 

_ i 

c3 

1 

_ i 

07  J 

The  prior  information  of  the  state  vector  after  two  time  steps,  q01  =  [q^^q,.]7 ,  is 
calculated  by  inserting  initial  information  from  stationary  states  as  follows 

^(q„:i  )  =  Q'  °- 

rs,(qrf,o:i)  o 

y  [  o 

Generalizing  Eq.  (0-104)-  Eq.  (0-106),  for  t  time  steps,  the  BFIM  from  prior  data, 
Jp( q0.,)  can  analytically  be  computed  by 


^(qto)  =  ^{-A^lnp(q0)} 

k=\ 

•%,  ®Ql  ®e;,‘+S  «  ®a 


k= 1 


-sl?l®Q?-s, 


t+ 2 


»^+cr+1 


0-91 


The  BFIM  elements  for  each  proposed  model  are  different  as  the  proposed  system 
models  vary  in  terms  of  state  variables  and  measurement  and  state  update  models.  In  this 


73/154 


regard,  the  following  subsections  explain  the  calculation  of  BFIM  elements  for  all  four  proposed 
models  in  more  detail.  The  BFIM  calculation  is  described  in  particular  for  the  fourth  system 
model  as  its  state  space  includes  both  discrete  and  continuous  variables. 


MODEL  1:  SYNCHRONIZED  ANS,  MN  SYNCHRONIZED  WITH  ANS,  KNOWN  SIGHT 
CONDITIONS 

Given  zikI]  N(ht  (q*.),crM  ,  the  Fisher  information  obtained  from  a  single  measurement 


from  the  i  -th  AN  is  obtained  by 

dh,( qkf  oh,(qk) 


0-92 


J'Z.i  ( fl/  ) 


dq *  dqk 


where  f  (q4)  =  -xk\  +  (mLv-yk)  ■  as  described  in  Eq.  (0-66).  The  gradient  of 
measurement/likeli^ood  PDF  is  calculated  by 


fyj q*). 

d<lk 

for  i  =  \, 2,... Nap 


ri,k 


'1x2  N„ 


i,k 


0-93 


and 


dfr(q*). 


0-94 


k v-^a) 

r,k 

A  i,x  k)  2(i-NAP)-\ 

L  cU2NFf  T 


2(i- 


l'lx2Aft.n  J 


foi  i  N ^ +1,N AP  +  2,...,N AP  +  N F 


where  rik  =  J(mix -xk)  +(miy-yk)  and  e\/m  is  the  1  x m  row  vector  whose  elements  are 


zero,  except  the  i  -th  element  which  is  equal  to  one.  Then,  the  measurement  Fisher  information 
matrix  from  a  single  measurement  available  at  time  step  t  from  the  i  -th  AN  is  calculated  by 
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0-95 


J z ,/  (q* ) _  cTz,, 


{mis~xk)2  {™^-xk){m,,y-yk) 

2  2 
ri,t  ri,k 

{m,,z~  Xk){m,,y-  yk)  [mi,y-yk) 


for  i  =  l,2,...NAP 

and 


Jz,i  (Qa  )  —  °"z,( 


ri,k 


'  i,k 


(mi,x-Xk)(mi,y-yk)  (mi,y~yk)~ 


’ i,k 

at 


’i,k 

bt 


A 

B 

D 


for  i  -  NAP  +  \,NAP  +  2,...,NAP  +  NFP 

where 

\2 


A  =  [- 


B  =  [- 


Mp).,  (mr}  -yk){m,,x-xk)  „2(/_JV„)l 


a 


'i,k 


K* -xk)(m.y-yk)  ^  ^  (m.y-yk) 

‘'IxlNn 


ri,k 


C=r(mi,x  Xk)  i,2(i-NAP)-l  {m:.v  -VA  )  i,2(i-*AP )  n 


i,k 

Nap)~ 

"N anx2N Fp  - 
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The  measurements  from  different  ANs  are  assumed  to  be  independent.  This  is 
reasonable  as  the  measurement  noise  comes  from  two  sources  which  are  the  receiver  noise 
and  the  multipath.  The  propagation  paths  from  the  ANs  to  the  MN  are  statistically  independent 
provided  that  the  ANs  are  sufficiently  separated  spatially  (which  may  be  violated  in  the  MIMO 
cases),  and  there  is  no  tunneling  propagation  or  keyhole  effect  at  the  receiver.  However,  in 
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some  cases,  the  multipath  occurs  close  to  the  MN  receiver;  hence,  if  two  ANs  are  at  similar 
bearings  relative  to  the  MN,  then  the  multipath  will  be  correlated.  However,  this  correlation  is  not 
considered  in  the  current  analysis.  Based  on  the  assumption  that  the  noise  in  the  receiver  is 
uncorrelated  as  the  synchronization  signals  from  the  ANs  have  different  spreading  codes; 
hence,  the  noise  subspaces  are  orthogonal.  Consequently,  the  BFIMs  from  additional 
statistically-independent  AN  signals  can  simply  be  added  as 


'  ’  AN  ‘'AN 


dhMk)'  dhMk) 
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1=1 


i= 1 


3q *  dqt 

To  comply  with  the  re-sorted  state  history  vector,  q()7  =  [qr/  07  qs.]r ,  the  elements  of 


Eq.  (0-102)  is  described  as  A’"  In  p(zk  |  qA)  is  the  upper  left  LdxLd  submatrix  of  Jz  (qA) , 
A’*  In p(zk  |  qk)  is  the  the  upper  right  LdxLs  submatrix  of  Jz  (qA ) ,  A^‘  In  p(zk  |  qA)  is  the 
lower  left  Ls  xLd  submatrix  of  Jz  (qA.) ,  and  A)]1  In p(zk  \  qk)  is  the  lower  left  Ls  xLs  submatrix 

of  Jz( q*)- 


To  calculate  the  BFIM  from  prior  information,  as  shown  in  Eq.  (0-107),  Q0  and  Qd  must 
be  defined  based  on  the  state  space  and  the  state  update  process  for  each  system  model.  In 
the  system  model  1,  the  vector  of  dynamic  states  qdk  is  only  described  by  MN  location 

variables  as  qdk=\xk  ykf  with  the  dimension  Ld  =  dim(qrf/)  =  2  .  Subsequently,  the  update 
process  for  dynamic  states  qrf0  and  qdl  is  described  by  q(/1-q,/0D  N(02xl,Qp ) ,  and  the  prior 
PDF  at  the  initial  point  is  defined  by 

0-99 

q0D  N(0 ,Qo) 

where  Q0  =  diag {cr^0  a;0  a^0  a^0  ...  a2mx^  (j  .  Q0  is  decomposed  into  the 

stationary  and  dynamic  portions  as 
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Qdo=dwg{(j;0  o-;0} 
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MODEL  2:  SYNCHRONIZED  ANS,  MN  NOT  SYNCHRONIZED  WITH  ANS,  KNOWN  SIGHT 
CONDITIONS 


For  this  system  model,  the  Fisher  information  from  the  measurement  of  i  -th  AN  is 


calculated  as 
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Similarly,  the  measurement  BFIM  from  all  MN-AN  links  is  calculated  by  the  sum  of  the 
FIM  from  each  MN-AN  link  as  stated  in  Eq.  (0-97).  It  must  be  noted  the  elements  of  FIM  from 
measurement  data  for  re-sorted  state  history  in  Eq.  (0-85)  are  obtained  in  a  similar  way  to  the 
procedure  explained  for  the  system  model  1 . 

To  calculate  prior  information,  the  initial  information  must  be  modified  by  adding  the 
initial  condition  of  the  range  offset,  given  by 


<\0UN(Q,Q0) 


where  Q0  =  drag  \  ax0  cr0  <0  cr 
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The  initial  condition  for  the  stationary  portion  remains  the  same;  however,  dynamic 
portion  is  modified  as  Qd0  =  diag\cr2x0  a20  cr2hQ  j .  The  dynamic  portion  of  state  vector  is 

defined  by  <\dk=  [ xk  yk  bkf  with  a  dimension  of  Ld  =3 .  The  update  information  is 
calculated  similar  to  Eq.  (0-97)-(0-107)  where  Qd  is  modified  for  newly  added  range  offset 
variables  as  Qd  =  diag\(Qp  ofj . 

MODEL  3:  UNSYNCHRONIZED  ANS,  MN  NOT  SYNCHRONIZED  WITH  ANS,  KNOWN 
SIGHT  CONDITIONS 

The  measurement  Fisher  information  is  calculated  similar  to  system  model  1,  as  stated 
in  Eq.  (0-91),  except  for  each  AN-MN  link,  a  single  range  offset  variable  is  added  to  the  state 

space,  and  h.(  qj  is  modified  by,  h,  (qt )  =  yj(mhx  -xkf  +  (miy  -  yk  )2  +  bik .  Therefore, 
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Consequently,  the  measurement  FIM  is  obtained  by 
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Similarly,  the  measurement  BFIM  from  all  MN-AN  links  is  calculated  by  the  sum  of  the 
FIM  from  each  MN-AN  link  as  stated  in  Eq.  (0-97).  It  must  be  noted  the  elements  of  FIM  from 
measurement  data  for  re-sorted  state  history  in  Eq.  (0-85)  are  obtained  in  a  similar  way  to  the 
procedure  explained  for  the  system  model  1 . 

The  range  offset  terms  make  the  rank  of  Jz  equal  to  one,  so  the  portion  of  Fisher 

information  that  relates  to  observation  is  not  invertible.  The  smooth  spatial  assumption  of  the 
MN  movement  is  necessary  to  make  the  FIM  invertible  and  the  SLAM  problem  observable.  To 
describe  the  prior  information  transformed  by  the  update  model,  the  initial  information  needs  to 
be  modified  by 

0-109 

q0UN(0,Q0) 


where  Q0  =  diag \ cr“0  cr ~  0  ...  ah  fi  crmy]  ()  ...  cr^  f.  The  initial 


condition  for  the  stationary  portion  remains  the  same  as  the  system  model;  whereas  the 
dynamic  portion  is  modified  by  Qd 0  =  diag^0  ay0  a^0  ...  cr2  0| .  In  this  case,  the 


dynamic  state  vector  qdk=\xk  yk  bJ7  has  a  dimension  of  Ld  =  2  +  Nan  .  The  update 
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information  is  calculated  similar  to  Eq.  (0-97)-(0-90)  where  Qd  is  modified  for  newly  added 


range  offset  variables  as  Qd  =  diag [Qp  Qb) . 


MODEL  4:  UNSYNCHRONIZED  ANS,  MN  NOT  SYNCHRONIZED  WITH  ANS,  UNKNOWN 
SIGHT  CONDITIONS 


As  mentioned  earlier,  one  approach  to  deal  with  mixed  NLOS/LOS  propagation  is  to 
introduce  two  different  noise  models  for  LOS  and  NLOS  propagation  conditions,  whereby  the 
transition  between  the  LOS  and  NLOS  states  is  modeled  with  a  binary  Markov  chain.  This 
modeling  turns  the  OWLS  to  a  jump  Markov  nonlinear  Gaussian  system  (JMNLGS).  Since  the 
FI  calculation  involves  taking  derivative  with  respect  to  all  state  variables,  it  cannot  handle 
discrete  values  parameters  such  as  binary  sight  conditions.  Therefore,  the  discrete  parameters 
need  to  be  marginalized  from  all  probability  densities. 

In  [liii],  [xii],  an  enumeration  method  is  proposed  to  approximate  the  desired  BCRLB  as 
the  expected  value  over  all  the  discrete  state  sequences,  referred  as  the  Enumer-BCRLB.  The 
method  initially  obtains  a  lower  bound  on  the  MSE  matrix  for  any  conditional  estimator  for  the 

continuous-valued  portion  of  q, ,  denoted  by  q' ,  based  on  the  calculation  of  the  conditional 
MSE  matrix,  given  by 
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MSE(qt  |  4)  □  Ep(q;_zu|4)  |[q;(z1:,  |  4)-q?][ttf  (z*  |  4)-q‘] 

^(q‘K) 1 

where  q‘(zl:(  |  s[.t)  denotes  the  conditional  estimator  for  a  particular  sight  state  sequence  s[.t 


and  /(q°  |  sj.,)  denotes  the  conditional  BFIM  for  q)  which  is  the  KxK  lower  right  submatrix  of 

J~ )  (qc0t  |  s[t ) .  Similar  to  the  procedure  explained  previously,  j(qct  |  can  be  calculated 

recursively  as  shown  in  Eq.  (0-75)  as  or  based  on  linearization  approximations  as  shown  in  Eq. 
(0-83).  Then,  to  evaluate  the  unconditional  MSE  lower  bound,  the  conditional  CRLB  is  averaged 
over  all  possible  discrete  state  sequences. 


81/154 


0-111 


Therefore,  the  total  BFIM  is  obtained  by  the  inverse  of  the  Enumer-CRLB.  It  can  be  seen 
that  calculation  complexity  of  Enumer-CRLB  algorithm  increases  tremendously  with  time, 
yielding  a  complexity  order  of  o{t. 2')  for  a  binary  sight  state  model  since  this  requires  the 
calculation  of  conditional  FI  Ms  for  all  possible  discrete  state  sequences.  In  addition  to 
computational  complexity,  it  is  stated  in  [xii],  [liv],  [Iv],  [[Ivi]  that  the  Enumer-BCRLB  can  be  an 
overly  optimistic  bound. 

To  obtain  a  tighter  bound,  [liv]-[lvi]  proposed  a  new  method  for  jump  Markov  linear 
Gaussian  systems  where  the  Enumer-BCRLB  is  modified  based  on  Monte  Carlo  (MC)  methods. 
Since  the  state  update  process  is  described  by  a  LJG  model  and  is  not  a  function  of  discrete 
states,  the  prior  information  for  continuous-valued  state  vector  can  be  calculated  analytically  as 


described  in  Eq.  (0-90).  However,  the  BFIM  from  measurement  data,  Jz{ q^)  is  numerically 
approximated  by  [Iv-lvi], 


0-112 


where  |q^|"  and  {zl7}"  are  the  independent  and  identically  distributed  vectors  generated  from 


p( q^,zl7)  for  N  Monte  Carlo  simulation  runs.  The  algorithm  to  compute  J z{(\{.t) 


approximately  is  adopted  from  [Ivi].  It  is  reminded  that  for  the  problem  at  hand,  q)  is  defined  by 
joint  vector  of  location  of  MN,  AN  locations  and  range  offsets  as  <\ct  =  [p,  b,  m]r . 

Unlike  the  Enumer-CRLB  ,  this  later  bound  does  not  require  explicit  calculation  of  all 
possible  discrete  state  sequences  which  reduce  the  complexity  to  the  order  of  o(2.N.t). 
However,  [Ivi]  shows  that  for  N  =  20000  ,  the  method  has  obtained  a  tighter  bound  than  Enumer- 
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CRLB  for  an  tracking  example;  for  an  nonlinear  SLAM-based  OWLS  with  the  high  dimensional 
state  space,  the  required  number  of  MC  runs  is  unknown  and  much  higher. 

Another  approach  for  calculating  BCRLB  for  jump  Markov  systems  is  suggested  in  [xii] 

that  is  based  on  deterministic  trajectory  of  discrete  states  sjt .  However,  it  can  be  overly 

optimistic  as  it  assumes  the  knowledge  of  discrete  state;  practically,  it  has  sufficed  requirements 
for  error  assessment  in  most  tracking  problems  [xii],  [xlv].  In  chapter  4,  it  will  be  shown  that  the 
BCRLB  obtained  based  on  the  knowledge  discrete  state  sequence  still  provides  a  reasonable 
bound  to  evaluate  the  result  of  the  proposed  algorithms  for  the  OWLS. 

Other  estimation  metrics  in  terms  of  the  BFIM 

An  accuracy  metric  is  important  in  that  it  gives  an  indication  of  confidence  in  OWLS 
solution  and  if  the  tracking  algorithm  works  successfully.  In  previous  section,  the  BFIM  is 
obtained  as  a  quantitative  metric  of  all  available  information  in  a  SLAM-based  OWLS.  The  BFIM 
comprises  all  information  regarding  the  parameters  of  interest  including  the  MN  trajectory  and 
nuisance  parameters  including  AN  locations,  range  offsets,  and/or  sight  conditions.  Since  the 
main  purpose  of  OWLS  is  MN  localization,  a  concise  accuracy  metric  which  only  focuses  on 
positioning  accuracy  can  be  more  useful  than  a  high  dimensional  BFIM.  In  this  regard,  Position 
Error  Bound  (PEB)  is  introduced  which  interprets  the  BFIM  in  terms  of  localization  accuracy. 
Moreover,  a  quantitative  confidence  metric  can  be  extracted  from  BFIM  known  as  confidence 
region  of  state  vector.  Not  only  in  terms  of  accuracy  and  confidence,  the  BFIM  introduces  an 
useful  tool,  known  as  Normalized  Estimation  Error  Squared  (NEES),  to  evaluate  OWLS 
solutions  in  terms  of  consistency.  In  following  subsection,  the  relation  of  these  three  important 
metrics  with  the  BFIM  is  explained  in  more  detail. 
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In  terms  of  localization  accuracy,  the  error  variance  of  the  MN  position  {xt,yt}  is  of 


interest,  which  is  the  upper  left  2by2  submatrix  of  Jtot  (q,) 


£{(p,-p,)(p,-p,f}^r*4  (q„r 


-12x2 


or 


E\(xt- xt )  +  (y>  - y, )“  P  tr\  Jtot  (q, ) 


where  tr 


J(qj  1  is  known  as  Position  Error  Bound  (PEB). 
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The  Normalized  Estimation  Error  Squared  (NEES)  for  the  state  q,  is  defined  as  [xlv] 

_  0-11 

=(q/-q,)  ^(qi)(q,-q() 

and  is  chi-square  distributed  ,  x~ ,  with  dim(q,)  degrees  of  freedom  under  the  hypothesis  that 

the  filter  is  consistent  and  approximately  linear-Gaussian.  In  an  estimation  problem,  particularly 
SLAM,  NEES  is  used  as  a  measure  of  filter  consistency  by  the  examination  of  average  NEES 
over  N  Monte  Carlo  runs  of  filters  [xliii],  [xlv];  [Ivii],  Given  N  simulation  runs,  the  average  NEES 
is  obtained  by 


N 


«  0-116 

2X 

i=\ 

The  quantity  of  Nsq  is  chi-squared  distributed  with  A^dim(q()  degrees  of  freedom  under 


consistent  filter  hypothesis.  The  95%  probability  concentration  region  for  sq  can  be  used  to  test 
whether  the  SLAM  filter  is  optimistic  if  the  error  rises  significantly  higher  than  bound,  or  is 
pessimistic,  if  it  tends  below  the  lower  bound.  Moreover,  can  be  asymptotically  estimated  by 

Afdim(q;) . 
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Confidence  Region  for  the  state  vector 

When  q,  is  the  unknown  parameter  to  be  estimated,  one  can  say  that  q,  must  be  within 
some  neighborhood  of  q,.  This  neighborhood  is  determined  by  the  confidence  region  of  the 
state  vector  q,  ,  as  defined  by  the  inside  of  “g-sigma”  ellipsoid, 

_  0-11\ 

=(q,-q,)  ^(q,)(q,-q,)  =  g2 

This  is  actually  a  hyper-ellipsoid  of  dimension  dim(q,)  where  the  semi-axes  are  g 
times  of  the  square  roots  of  the  Eigen  values  of  [j(q,)]  ' .  For  an  efficient  estimator,  the 
normalized  error  sq  must  follow 
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where  Q  represents  the  small  tail  probabilities.  For  a  one-sided  probability  region  of  95%,  the 
ellipsoid  boundary  will  be  used  to  assess  the  SLAM  algorithm  error  performance  in  comparison 
with  the  efficient  estimator. 

Mutual  information  (MI),/(X,7),  measures  the  mutual  dependence  of  two  random 
variables,  X  and  7.  If  X  and  7  are  the  input  and  output  of  a  channel,  respectively,  I(X,Y ) 

does  not  make  any  assumptions  regarding  the  receiver  processing  of  7,  but  requires  the 
source  entropy.  It  can  be  shown  that  if  the  signal  to  noise  of  each  degree  of  freedom  (DOF)  is 
small,  then  I(X,Y)  is  independent  of  the  source  distribution  [Iviii],  [lix],  [lx].  Hence,  an 
asymptotic  assessment  of  the  tracking  performance  based  on  the  mutual  information  is  possible 
for  any  source  distribution.  We  can  select  a  source  distribution  of  X  that  facilitates  the 
mathematics  and  results  in  the  most  convenient  expressions.  Assuming  a  low  SNR  per  signal 
DOF,  such  that  the  source  distribution  becomes  asymptotically  irrelevant,  then  it  is  possible  to 
consider  the  symmetry  of  the  problem  such  that  the  mutual  information  can  be  assessed  for  a 
specific  value  of  X .  However,  instead  of  determining  the  overall  I(X,Y ) ,  where  the  probability 
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structure  of  X  needs  to  be  specified,  we  can  consider  the  simpler  problem  of  determining 
I(X,Y\X  =  A)  or  just  I(Y,A ) .  Mutual  information  quantifies  how  much  information  of  a  specific 
value  of  A  is  contained  in  Y .  In  other  words,  determine  the  statistics  of  the  optimum  estimator 
of  A  based  on  the  data  in  Y  which  is  expressed  as  A(Y ) .  The  quantity  of  information 
regarding  the  parameter  A  that  is  contained  in  Y  is  also  given  by  the  Fisher  information 
denoted  in  this  case  as  J(A;Y ) .  Hence,  the  quantity  that  is  used  to  assess  the  estimation  and 
tracking  performance  is  the  effective  signal  to  noise  given  as  A2J(A;Y) . 

Guo  et  al  [xi]  revealed  an  interesting  relation  between  the  mutual  information  and  the 
MMSE.  Guo  theorem  says  that  regardless  of  the  input  statistics  the  relation  between  mutual 
information  and  MMSE  in  an  additive  Gaussian  channel  can  be  described  by 

dl  (X,  Y;SNR)  1 

v  - '  =  -  MMSE 

dSNR  2 

where  SNR  is  the  signal-to-noise  ratio  of  the  channel;  I (X;7)  is  the  input-output  mutual 

information  in  nats,  and  the  MMSE  is  the  minimum  mean  square  error  of  input  estimation  given 
output  as  a  function  of  SNR.  Eq.  0-118)  reflects  a  strong  relevance  of  mutual  information  to 
estimation  and  filtering  which  can  provide  a  non-coding  operational  application.  The  input-output 
mutual  information  demonstrates  how  much  information  can  be  transferred  through  a  channel 
reliably  given  a  certain  input  signaling,  while  the  MMSE  quantifies  how  efficiently  each  input 
sample  can  be  recovered  using  the  channel  output.  The  significance  of  this  relation  is  that  this  is 
not  only  simple  and  intriguing  but  is  also  independent  of  input  distribution  which  holds  under 
arbitrary  signaling  conditions  and  the  broadest  setting  of  the  Gaussian  channel,  including  both 
discrete-time  and  continuous-time  channels  in  both  scalar  and  vector  versions.  Two  explanatory 
proofs  for  Guo’s  theorem  are  presented  in  Appendix  A  where  the  first  proof  is  inspired  by 
independence  of  Ml  from  input  distribution  in  low  SNRs  and  incremental  channel  concept  and 
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the  second  proof  uses  the  relation  between  FI  and  differential  entropy  in  de  Bruijn’s  identity  [xi], 

[Ixii]. 

The  intuitive  Guo’s  theorem  introduces  a  novel  application  of  another  information  metric, 
that  is  mutual  information,  in  tracking  and  estimation  problems  [Ixi].  However,  the  Guo  theorem 
only  proposed  a  Ml  application  in  a  linear  Gaussian  problem  for  random  stationary  parameter.  In 
this  section,  it  is  explained  how  the  mutual  information  concept  can  be  applied  for  a  SLAM- 
based  OWLS  which  includes  not  only  random  stationary  parameters  such  as  FP  locations  but 
also  random  dynamic  variables  such  as  the  MN  location.  The  elements  of  Ml,  their  contribution 
to  stationary  and  dynamic  parameter  estimation  are  depicted  via  a  Venn  diagram  based  on  a 
Markovian  structure  of  OWLS. 

We  define  a  SLAM-based  wireless  localization  problem  where  the  MN  moves  within  an 
unregistered  network  and  receives  wireless  observables  from  a  mixture  of  stationary  ANs  with 
unknown  or  known  positions.  The  MN  trajectory  is  assumed  to  be  random  but  smooth  in  the 
way  that  the  motion  dynamical  process  is  modeled  by  a  Markovian  process.  At  any  discrete  time 
instant  /eD  ,  the  state  variable  q,  includes  stationary  parameters  such  as  FP  locations  m  , 

and  dynamic  variables,  such  as  the  MN  location,  and  other  nuisance  parameters  such  as  range 
offset  qdt  =[p,  b,]r,  is  denoted  as 


0-120 


Mutual  information,  /( q,,q,)  is  a  measure  of  information  which  the  sufficient  statistic  or 
estimator  ,q; ,  can  provide  for  estimating  q;  in  terms  of  shared  entropy  [Ixii].  The  first  element  of 
Ml  is  defined  by  Markovian  process  between  subsequent  dynamic  states.  Note  that  the 
stationary  portion  of  state  vector  m//;  can  be  also  considered  as  a  dynamic  variable  with  a  zero 

mean  and  variance  Markov  process.  Figure  0-20  shows  a  Venn  diagram  of  entropy  change  in 
two  and  three-step  Markov  process.  The  entropy  of  states  are  shown  by  color-filled  ellipse  areas 
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and  their  common  information  or  Ml  shown  by  overlapped  areas  between  entropies.  For 
example,  the  green  left  area  in  left-side  depicts  the  maximum-decreased  entropy  of  qdt  if  the 
previous  state  is  known  as  a  result  of  filtering  in  a  online  processing,  while  green  area  in  right 
side  diagram  shown  the  maximum-decreased  entropy  of  qdt  if  previous  and  next  states  are 
known  as  a  result  of  smoothing  in  a  batch  processing. 


Since  the  underlying  Markov  process  (q^M  — »  qdJ  — »  qdt+l )  evolves  unaware  of  the 
observations,  when  at  the  time  instant  t  ,  the  observation  z,  depends  only  on  the  state  at  time 
instant  t  .  This  means,  given  qdt ,  the  observation  at  the  time  instant  t  ,  z,and  q^are 
independent  so  that  qdt ,  z,and  qdt_x  form  a  Markov  Chain  as  qdt_x  — »  qd  T  — >  zf  or 

Qdr-i  that  is  equivalent  with  I(qdt_1;zt  \  qdt)  =  0  .  Note  that  this  independency  is  also 

applied  to  the  whole  state  vector  q,.  ,  since  the  stationary  portion  of  state  vector  can  also  be 

treated  as  dynamic  variable  as  earlier  discussed.  Therefore,  the  second  element  of  Ml  is 
defined  by  measurement  data;  however,  this  independency  property  must  be  considered  when 
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combining  the  Ml  elements.  Figure  0-21  shows  entropy  change  in  presence  of  measurement 
data  and  Ml  dependency  in  this  Markov  process. 


Figure  0-21  Venn  diagram  of  the  entropy  change  in  presence  of  measurement  data 

With  Markovian  chain  assumption  between  states  and  also  state-measurement,  two 
important  points  is  taken  into  account  in  plotting  Venn  diagrams.  First  /( qM;z,  |  q, )  =  0  ,  then 

0-121 

/z(qM;z,)  =  /z(qM;z,|q,) 

where  h(.)  denotes  the  entropy.  It  means  that  not  only  /(qM;z/)</(qM;q<)  but  also 

(qMnz,)c(qMnq,). 

Second,  since  p{(\,  |  qM,zM)  =  p( q,  |  qM) ,  it  is  inferred  that  given  the  previous  states, 

previous  measurement  doesn’t  add  new  information  for  the  current  and  future  states,  which 
complies  with  Markovian  assumption. 
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Therefore,  the  maximum  mutual  information  for  state  estimate  at  time  step  t  in  terms  of 
nats  or  bits  is  obtained  by  combining  the  Ml  elements  with  considering  the  Markov  properties  , 
given  by 

J(q,  »q,-i)  ^  J(q,  ;q,-p  z,)  =  /(qr  ;z,)  +  /(q,  ;q^i)-/(q(_i;z()  0-122 

=  J(q,;z,)+/(q,;q,_i|z,) 


Figure  0-22  Venn  diagram  of  the  Ml  and  entropy  change  in  a  Bayesian  problem 

Finally,  the  entropy  reduction  and  total  Ml  increase  procedure  of  three  subsequent  states 
in  a  SLAM  problem  is  illustrated  in  Figure  0-22  which  fulfills  these  two  criteria.  In  an  online 

SLAM  problem,  Eq.  (0-121)  determines  the  upper  bound  of  information  for  q, ;  while  in  the  full 
SLAM  problem,  the  added  information  from  future  states  must  be  considered,  as  shown  Figure 
0-20,  where  the  decreased  entropy  of  q,  after  taking  into  account  the  information  from  is  q,+1 

and  zt+l  illustrated  by  the  pure  green  area.  However,  the  Ml  calculation  for  state  variables  given 

in  Eq.  0-121)  can  be  difficult  in  comparison  to  MMSE  calculation;  it  gives  a  global  criterion  for 
any  state  statistics  and  channel  noise  distributions  in  units  of  bits  or  nats.  The  more  detailed 
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pros  and  cons  of  Ml  application  in  estimation  problems  will  be  discussed  in  next  subsection  by 
some  examples. 

Guo’s  theorem  applications  and  its  limitation  for  the  SLAM-based 
OWLS 

The  main  advantages  of  Ml  application  to  estimation  problems  in  comparison  with  FI  or 
Bayesian  CRLB  is  that  it  provides  a  global  measure  of  shared  information  between  the  state 
vector  and  observables  and  prior  information  regardless  of  linearity  or  Gaussianity  of  the 
problem.  The  Guo’s  theorem  has  built  the  bridge  between  estimation  theory  in  terms  of  MMSE 
and  information  theory  in  terms  of  Ml  under  certain  conditions.  However,  Guo’s  theorem  and  Ml 
calculation  have  their  own  disadvantages  which  lead  us  to  mainly  focus  on  the  FIM  for  our 
SLAM-based  wireless  localization  problem.  Following  example  are  designed  to  illustrate  the  Ml 
and  Guo’s  theorem  limitations  for  OWLS  application  in  more  detail. 

Example  01:  Linear  Gaussian  random  input  in  additive  Gaussian  channel 

Consider  a  simple  channel  model  with  continuous-valued  input  and  noise  process 

r  0-123 

y - ySx  +  n 

where  n,xU  N( 0,1) .  Therefore,  the  output  distribution  is  given  by 

0-124 

yU  iV(0,£  +  l) 

For  a  Gaussian  channel  with  linear  Gaussian  input,  the  input-output  mutual  information 
in  terms  of  nats  is  derived  as 

/(x;  y)  =  H(y)~  H(y  \  x)  =  j  ln(2xe(d  + 1))  -  ^  ln(2^e)  =  ^ ln(£  + 1) 

Figure  0-23  shows  the  mutual  information  calculated  numerically  from  entropy  and 
analytical  result  in  Eq.  (0-124).  Although  the  Ml  calculation  from  entropy  is  trivial  when  the 
closed  form  is  obtainable,  both  plots  are  presented  in  order  to  observe  possible  errors  of 
numerical  integration.  From  Guo’s  theorem,  it  is  inferred  that  the  mutual  information  of  a 
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Gaussian  channel  is  concave  function  of  SNR  as  shown  in  Figure  0-21 .  Moreover,  the  mutual 
information  can  be  bounded  as 


mmse(SNR)  <  ^  I(SNR)  <  mmse{ 0)  =  var(X) 


Ml  vs  SNR  for  linear  gaussian  input 


0-126 


SNR(dB) 


a.  linear  scale  SNR  b.  SNR  in  dB 

Figure  0-23  Mutual  information  versus  SNR  in  a  linear  Gaussian  channel  with  a  Gaussian 
input 


The  MMSE  of  the  input  estimate  given  the  output  is  computed  as  the  averaged  mean  of 
the  posterior  PDF  : 

fY\x(y\x)fx(x)  °'127 


fx\y(X\y)  =  - 
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Since  the  variance  of  posterior  PDF  is  constant  for  all  value  of  y  ,  the  MMSE  is  given  by 

1  0-128 


1 +  S 
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GUO’s  theorem  is  confirmed  not  only  by  analytical  results  in  Eq.  0-128)  for  both 
continuous  Gaussian  signal  and  noise  in  a  linear  system  but  also  by  numerical  calculation  as 
shown  in  Figure  0-22: 


S"*A=X-mmse(x\ S,y)  =  —  °'129 

dS  2  \  +  S 

Small  difference  between  analytical  and  numerical  results  for  SNR<10dB  in  Figure  0-24 
is  mainly  due  to  numerical  error  in  Ml  first  derivative  calculation.  However,  it  is  not  mentioned  in 
[xi],  Eq.  (0-128)  is  only  valid  when  the  variance  of  x  is  normalized  to  1. 


MMSE  vs  SNR  for  linear  gaussian  input 


Figure  0-24  MMSE  and  Ml  versus  SNR  in  linear  Gaussian  channel  with  Gaussian  input 


Example  02  :  Nonlinear  discrete  random  input  in  an  additive  Gaussian  channel 

Range  and  TOA  measurements  provide  observables  for  tracking  algorithms  that  have 
nonlinear  relation  with  state  space.  In  this  case,  the  verification  of  Guo’s  theorem  for  nonlinear 
system  is  investigated.  Even  if  the  input  x  is  Gaussian  ,  the  output  PDF  through  nonlinear 
function  is  not  Gaussian.  One  of  main  disadvantages  of  Ml  application  in  tracking  problems  will 
be  the  difficulty  in  generating  the  PDF  for  a  general  nonlinear  process  with  some  correlation.  In 
this  regard,  numerical  calculation  is  generally  solution  such  as  approximation  of  the  continuous 
PDF  of  input  by  a  discrete  PDF.  Consider  a  input-output  nonlinear  system  with  additive 
Gaussian  noise  channel: 
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0-130 


y  -  \[S  f(x)  +  n 

M  M 

where  nU  jV(0,1)  and  xU  Y,P‘s{x~x) as  I>  =  '  ■ 

i= 1  i=l 

In  this  example,  a  Gaussian  PDF  of  input  x  is  approximated  by  a  discretized  version  of 
Gaussian  PDF  with  cell  width  A.  To  apply  Guo’s  theorem,  it  is  also  required  to  normalized  the 

variance  of  /(x) ,  (Jfx ;  therefore  the  total  SNR  is  stated  by  \[8cj ^  where 


v  =  \[S<j 
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+  n 
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the  output  distribution  is  given  by 
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The  posterior  PDF  is  calculated  as 

fY\X(y  I  x)fgM(g(x)) 
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and  the  conditional  mean  estimate  is  given  by 

£0)  =  J  gO  )/gWLv  (gO)  I £*) 
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The  averaged  variance  of  posterior  PDF  which  is  equal  to  MMSE  is  given  by 
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MMSE(g(x )  |  y)  =  J  fY  (v)var(  /;,(t)|>  (g(x)  |  y)dy 
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SLAM-Based  Opportunistic  Wireless  Localization 

OWLS  is  proposed  as  a  solution  for  localization  in  environments  where  GPS  fails  or  the 
local  wireless  network  signaling  is  not  intended  for  positioning  purposes.  It  is  shown  that  the 
SLAM  can  introduce  a  systematic  approach  of  accounting  for  the  various  disparate  sources  of 
information  to  track  MN  trajectory  in  an  unknown  environments.  The  significance  of  SLAM 
application  in  an  OWLS  problem  is  that  it  can  incorporate  all  available  information  from  the 
motion  process,  observables  and  prior  knowledge  to  track  the  MN  trajectory  while  jointly 
estimates  the  unknown  environment  map  parameters  including  FP  locations  and  other  nuisance 
parameters  such  as  range  offset  due  to  unsynchronized  reception  and  multipath  effects. 

Chapter  2  introduced  the  OWLS  system  parameters  which  based  on,  Chapter  3  presented  four 
system  models  for  an  OWLS  according  to  range  measurements  data.  It  is  shown  in  Chapter  2 
that  optimal  recursive  Bayesian  solution  for  an  OWLS  requires  the  complete  posterior  PDF  of 
the  system  state  to  be  determined  as  a  function  of  time.  However,  a  closed  form  solution  for  an 
OWLS  cannot  be  formulated  due  to  system  nonlinearity.  Despite  the  absence  of  closed  form 
analytical  solution,  Bayesian  FIM  was  formulated  for  all  four  proposed  system  models  as  the 
best  achievable  performance  of  the  OWLS  in  Chapter  3. 

This  chapter  presents  two  suboptimal  Bayesian  solutions  for  an  OWLS  based  on  EKF 
and  PF,  known  as  EKF-SLAM  and  FastSLAM.  The  proposed  suboptimal  EKF-based  and  PF- 
based  solutions  are  implemented  for  four  proposed  OWLS  models  ranging  from  a  synchronized 
scenario  with  known  sight  conditions,  as  in  system  model  1,  to  where  there  is  no  knowledge 
about  the  MN  trajectory,  AN  locations  in  a  non-stationary  multipath  situation.  Figure  0-25 
illustrates  a  flow  chart  of  how  proposed  algorithms  interacts  with  system  models.  This  step-by- 
step  modeling  provides  a  deep  understanding  of  how  much  information  each  additional  variable 
contributes  to  solution.  BFIM  and  its  derivatives’  bounds  (BCRLB,  NEES,  PEB  and  confidence 
region),  formulated  in  Chapter  3,  are  used  as  benchmarks  for  comparison  of  implemented 
suboptimal  algorithms  and  assessment  of  the  effect  of  introduced  system  parameters. 
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Performance  bounds  are  also  used  to  investigate  the  effect  of  introduced  approximation  in 
suboptimal  solutions  in  terms  of  consistency  and  efficiency. 


Joint  MN  location  &  AN 
locations  state  vector: 
represented  by  a 
Gaussian  PDF 


MN  location: 
represented  by  particles 
AN  locations: 
represented  by  a 


Joint  MN  location,  range 
offset,  &  AN  locations 
state  vector:  represented 
by  a  Gaussian  PDF 


Joint  MN  location  and 
range  offset  state  vector: 
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Joint  MN  location,  range 
offset,  &  AN  locations 
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MN  location: 
represented  by  particles 
Joint  AN  locations  and 
range  offsets  state 
vector:  represented  by  a 
Gaussian  PDF 


Same  as  system  model 
3,  LML  as  sight  state 
detector 


Same  as  system  model 
3,  LML  as  sight  state 
detector 


Figure  0-25  Flowchart  of  OWLS  algorithms  according  to  system  models 

based  and  PF-based  solutions  for  an  arbitrary  trajectory  of  a  single  MN  in  a  rectangular  room. 
Figure  0-26  illustrates  the  Matlab-based  user  interface  to  simulate  an  OWLS  scenario.  Initially, 
an  arbitrary  MN  trajectory  is  generated  by  choosing  the  spline  points  shown  by  yellow  squares 
in  Figure  0-26.  The  MN  trajectory  is  depicted  by  blue  line  and  the  MN  locations  where  it  receives 
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new  measurements  and  updates  its  position  is  shown  by  red  dots.  The  number  and  the  location 
of  FPs  (shown  by  a  red  square  box)  and  APs  (shown  by  a  green  square  box)  can  be  chosen  by 
user  visually  over  the  trajectory  plot.  To  generate  both  NLOS  and  LOS  conditions  along  the 
simulated  trajectory,  a  wall  is  simulated  which  can  be  created  and  moved  by  user,  depicted 
green  rectangular  in  Figure  0-26 


>> 

\ 

\ 

\ 

01 - - - - - - - - - - 

01  23456789  10 

X  (meters) 


Figure  0-26  User  interface  to  simulate  MN  trajectory  and  locations  of  FPs  and  APs 


In  next  stage,  the  toolbox  generates  measurements  according  to  the  selected  system 
model  and  simulates  EKF-based  and  PF-based  solutions  simultaneously  as  illustrated  in  Figure 
0-25.  The  95%  confidence  regions  obtained  by  approximated  posterior  PDF  by  tracking 
algorithm  and  BFIM  formulation  from  Chapter  3  are  generated;  where  the  former  is  shown  by 
ellipsoids  with  solid  line  (red  for  FP  location  estimate  and  blue  for  the  MN  location  estimate), 
and  the  latter  is  shown  by  ellipses  with  dashed  line.  The  figure  also  compares  the  true  MN 
trajectory  with  the  estimated  MN  trajectory,  shown  by  black  thick  line  with  yellow  dots 
representing  the  MN  location  estimates,  and  FP  location  estimates  trajectory  are  illustrated  by 
green  dots.  The  MN  localization  errors  can  be  appreciated  by  looking  at  the  black  lines  that 
connect  the  true  and  the  estimated  MN  positions.  The  final  estimates  of  MN  location  and  FP 
locations  are  shown  by  purple  donuts  where  the  initial  estimates  are  shown  by  black  stars.  The 
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start  point  of  the  MN  trajectory  can  be  detected  in  EKF-based  results  as  the  red  point  which  is 
centered  by  two  blue  solid  ellipsoids  where  the  larger  one  represents  the  95%  confidence  region 
obtained  from  initial  conditions.  A  listing  of  figure  legend  descriptions  is  summarized  in  Table 


0-2. 


0  2  4  6  8  10 
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a) 
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b) 


Figure  0-27  Example  of  EKFSLAM/FastSLAM  tracking  results 
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Table  0-2  Figures’  legend  guideline 


Legend 

Description 

Blue  line 

The  MN  true  trajectory 

Red  dots  on  blue  line 

MN  estimate  update 

Yellow  square 

Spline  points  for  simulating  the  MN  trajectory 

Green  box 

AP  location 

Red  box 

FP  location 

Green/blue  rectangular 

Wall 

Black  star 

Initial  estimate 

Purple  donut 

Final  estimate 

Green  dot 

FP  estimate 

Yellow  dots 

MN  location  estimates 

Red  dots 

Particle-conditioned  FP  estimates 

Black  dots 

MN  location  particles 

Red  ellipse  with  solid  line 

95%  confidence  region  from  posterior  PDF  for  the  FP 
estimates 

Blue  ellipse  with  solid  line 

95%  confidence  region  from  posterior  PDF  for  the  MN 
estimates 

Red  ellipse  with  dashed  line 

95%  confidence  region  from  BFIM  for  the  FP  estimates 

Blue  ellipse  with  dashed 
line 

95%  confidence  region  from  BFIM  for  the  MN  estimates 

System  model  1:  synchronized  ANs,  MN  synchronized  with  ANs,  known  sight  conditions 

In  the  first  case,  the  MN  receives  signal  from  an  arbitrary  number  of  stationary  and 
synchronized  ANs  including  both  FPs  and  APs.  The  MN  clock  drift  from  ANs  clock  is  not  an 
issue  since  the  MN  is  synchronized  with  ANs  initially  in  calibration  stage  and  then  via  a  global 
reference  clock  such  as  GPS  clock  as  the  MN  moves.  As  earlier  described  in  Section  0  ,  the 

location  state  MN  p,  =  \x,,  ytf  is  modeled  as  a  first-order  Markov  process  ruled  by 

P,=Pm+v„  where  the  vf Q ,N(Q,Qp) .  Covariance  matrix  of  motion  process,^,  can  be 
chosen  commensurate  with  the  velocity  variance  of  the  mobile  user.  However,  a  more 
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informative  motion  model  can  be  replaced  from  output  CV  or  I  MU  sensors,  whenever  available. 
The  current  state  vector  of  unknown  variables  is  denoted  by  qr  =  [pr  m/p]  ,  where  m  fp  is  the 

stacked  vector  of  the  FP  locations,  defined  by  \mN  m,,  mM  mN  Y.  The 

L  jyAP+ 1 » x  ^AP+X  >y  yV AN  AN  >y  J 

estimate  of  q,  must  be  extracted  from  the  history  of  measurements  by  the  time  step  t , 

z1:,  □  {zj  z2  ...  zf} .  zk  =  jzj k  z2k  ...  zNm  k\, k  =  ,  denotes  range  measurements 

at  the  time  step  k  from  all  ANs  in  line  of  sight  with  the  MN.  It  is  assumed  that  the  ANs  provides 
LOS  measurements  or  if  not,  the  sight  condition  are  known  to  the  MN  and  NLOS  measurements 
are  discarded  from  the  measurement  set.  The  suboptimal  estimate  of  q,  is  obtained  by 

evaluation  of  the  posterior  PDF,  p(qt  |zw),  using  EKFSLAM  and  FastSLAM. 


EKFSLAM  solution  for  case  1 

Based  on  system  model  1  outlined  in  Section  0,  EKFSLAM  for  this  case  approximates 
the  posterior  PDF,  p(qt  |  zH)  =  p( p;,m/(,  |  zw) ,  by  a  joint  Gaussian  PDF  whose  the  mean  vector 
and  the  covariance  matrix  are  given  by 


q,i,= 


P/|r 

111 


fp,' 


P/|r-l 
AJPJ- 1 


+ 


K,[z,-h(p, 


0-136 


+K,D,K,r 

where  q,M  =  [p,M  t  ,  J  and  Z,M  are  the  mean  and  covariance  matrix  of  apriori  PDF, 
stated  by 


Pi|r~l 

i 

T 

7 

_ i 

q#-l  = 

Vl 

= 

1 

3- 

7 

i _ 

i 

3> 

7 

i _ 

1^1  -  Emm  +FtQp,,F 

and 
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0-137 


F  = 


1  0  0  0  •••  0 

0  1  0  0  •••  0 


2+2  Npp 


D  t=HZt[t_lHT  +  Rt 

where  H  is  the  Jacobian  of  range  measurement  function,  evaluated  at  q,^  =[p/;_1,m(_l]7  , 


previously  defined  in  Eq.  (0-92)-(0-93). 


FastSLAM  solution  for  case  1 

The  FastSLAM  solution  is  based  on  the  Rao-Blackwellization  theorem  where  the  joint 
posterior  PDF  can  be  factored  into  the  MN  location  components  and  a  conditional  ANs  location 
component  as  described  in  Eq.  (0-27).  It  approximates  the  posterior  PDF  by  set  of  particles 
where  each  particle  of  posterior  PDF  represents  the  history  of  MN  locations  which  additionally 
contains  the  Gaussian  parameters  (i.e.  the  mean  and  covariance  matrix)  of  FP  locations 
estimates.  The  FastSLAM  for  system  model  1  adopts  the  same  procedure  explained  previously, 

Eq.  (0-27)-(0-31)  where  the  measurement  function  /?,(.,.)  is  described  by  the  range 

measurement  function  stated  in  Eq.  (0-66). 

The  performance  of  the  SLAM  algorithms  is  assessed  by  a  simulated  localization 
scenario  shown  in  Figure  0-28.  The  MN,  while  moving  within  a  environment  of  10m  by  10m  , 
localizes  itself  and  maps  FPs  using  the  range  measurement  obtained  from  all  ANs  in  the  LOS 

condition.  A  smooth  MN  trajectory  jp,  j1  is  simulated  within  the  allowed  2-D  space  (walls  and 

the  AN  locations  are  not  permitted).  Sight  conditions  are  generated  by  accessing  the 
LOS/NLOS  map,  as  shown  in  Figure  0-28  for  four  ANs,  according  to  the  simulated  MN 
trajectory.  It  can  be  noticed  that  most  of  the  area  is  always  covered  by  at  least  2  ANs  which  their 
locations  are  not  necessary  known  to  the  MN.  This  makes  the  localization  task  particularly 
complicated  and  requires  MN  tracking  in  order  to  avoid  large  errors  due  to  poor  network 
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cooperation.  Nevertheless,  it  will  be  shown  that  the  Bayesian-based  SLAM  algorithms  will 
provide  good  localization  performances. 
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Figure  0-28  Sight  condition  map  for  4  ANs  along  with  the  MN  trajectory  in  a  2D  10m  by 
10m  room:  NLOS  propagation  area  for  every  AN  is  illustrated  by  gray  color 


Measurements  are  simulated  according  to  Eq.  (0-66)  and  EKF-SLAM  and  FastSLAM 
algorithm  are  implemented  for  different  scenarios:  4  APs,  3  FPs-1  AP,  and  4  FPs  to  investigate 
the  effect  of  knowledge  of  the  AN  location(s)  over  localization  accuracy. 

A  single  Monte  Carlo  simulation  is  carried  out  for  the  tracking  example  illustrated  in 
Figure  0-29  based  on  system  model  1,  referred  as  case  1.  For  better  visual  purposes,  the  MN 
path  has  been  generated  smoother  and  shorter  for  only  100  points  and  the  two  sided  95% 
confidence  region  ellipse  for  the  MN  location  estimate  is  plotted  for  every  10  time  step.  This 
figure  compares  the  true  trajectory  (blue  line  with  red  square  representing  the  sampled  actual 
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MN  location)  with  the  estimated  one  (black  thick  line  with  yellow  dots  representing  the  MN 
location  estimates)  obtained  by  EKF-SLAM  in  left  side  and  FastSLAM  in  right  side.  The  MN 
localization  errors  can  be  appreciated  by  looking  at  the  lines  that  connect  the  true  and  the 
estimated  MN  positions  in  Figure  0-29(a)-(e).  It  must  be  noted  in  all  simulations  the  true  MN 
location  is  inserted  as  initial  condition  of  MN  trajectory  which  can  fix  the  origin  of  coordinate 
system;  however,  the  coordinate  system  still  has  rotational  DOF  unless  another  known-located 
point  fixes  the  coordinate  system  such  an  AP  or  a  FP  can  be  located  with  good  approximation 
before  the  algorithm  diverges.  The  number  of  ANs  in  NLOS  conditions  are  plotted  in  Figure  0-30 

for  each  position  of  the  MN  along  the  path.  It  can  be  seen  that  for  time  interval  [22,3 1] ,  half  of 

ANs  are  blocked  by  the  wall  and  FP  locations  estimates  still  have  high  uncertainty,  leading  to 
large  MN  localization  errors. 


105/154 


Y  (meters)  ^  (meters)  Y  (meters) 


EKFSLAM  with  0  Anchor  Points, 4  Feature  Points 


FastSLAM  with  0  Anchor  Points, 4  Feature  Points 


X  (meters) 

c) 


0  2  4  6  8  10 


X  (meters) 

e) 


X  (meters) 

d) 

FastSLAM  with  4  Anchor  Points, 0  Feature  Points 


'X 

4 ^ 

)  TAT 

Wl  ( 

/A 

i 

Y  ^ 

/  u  ) 

s' 

0  2  4  6  8  10 

X  (meters) 


f) 


Figure  0-29  Example  of  MN  tracking  and  ANs  mapping  for  system  model  1 
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a)  4  FPs  b)  3  FPs-  1  AP 

Figure  0-30  Number  of  ANs  (APs/FPs)  in  NLOS  condition 


The  95%  confidence  region  for  the  MN  location  and  FP  locations  calculated  based  on 
corresponding  submatrixes  of  [Jto,  (qw)]  '  to  their  state  variables.  The  semiaxes  of  the 
confidence  region  ellipse  is  g  times  of  square  roots  of  the  Eigen  values  of  this  corresponding 
submatrix.  For  example,  the  ellipse  formula  for  the  MN  location  estimate  is  calculated  as 


(p,-p,)r-4  (p„)(p,-p,)  =  r 


0-138 


where  Jtot{ p,)  = 


MN 

J2x2 


and 


jtoX 


MN 


2x2 


is  the  corresponding  submatrix  of 


[ J(o;  (q|7 )]  'to  the  MN  location  variables. 


For  a  single  run  MC  and  2D  location  variable 


g  -  F/2( 0.95)  =  2.447  ,  where  F  2‘  is  the  inverse  function  of  second  order  chi-square  CDF.  The 

confidence  region  obtained  by  the  BFIM  is  also  compared  with  the  confidence  region  provided 
by  the  covariance  matrix  of  the  approximated  posterior  PDF  from  EKFSLAM  or  FastSLAM.  For 
FastSLAM,  the  posterior  PDF  covariance  matrix  of  the  MN  location  estimate  is  approximated 
from  distribution  of  particles,  and  for  the  AN  location  estimate,  is  approximated  by  spread  of  Np 
FP  estimates  obtained  in  each  EKF  run. 
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For  a  single  MC  run,  it  can  be  seen  that  FastSLAM  and  EKFSLAM  performance  in  terms 
of  MN/FP  location  estimates  are  comparable  for4-FP  scenario.  However,  the  concentrated 
estimates  of  FP  locations  in  red  dots  illustrates  that  FastSLAM  becomes  overconfident.  This  is 
also  observed  from  the  large  difference  between  the  confidence  region  obtained  from  BFIM  and 
estimated  FP  covariance  matrix.  This  is  due  to  inability  of  FastSLAM  to  maintain  particle 
diversity  over  long  periods  of  time.  The  effect  of  single  AN  location  knowledge  over  EKFSLAM 
and  FastSLAM  performance  can  be  observed  by  comparing  Figure  0-29  (a)-(b)  to  Figure  0-29 
(c)-(d).  Finally,  in  Figure  0-29  (e)-(f)  where  it  is  assumed  that  there  is  no  ambiguity  about  FP 
location,  the  EKF  performed  far  superior  regarding  uncertainty  estimation;  however,  in 
FastSLAM,  uncertainty  estimation  degrades  as  a  function  of  trajectory  length. 

For  all  ANs,  Figure  0-31  shows  the  actual  MN-AN  distance  and  the  corresponding 
estimates  by  EKFSLAM  and  FastSLAM.  For  the  NLOS  area,  the  distance  is  estimated  based  on 
the  last  FP  location  estimate  and  current  MN  location  estimate.  It  can  be  seen  that  however  the 
MN  localization  error  can  reach  to  2  meters,  the  ranging  error  remains  lower  than  1  meter, 
except  NLOS  areas  of  AN  03  and  AN  04  where  the  FP  is  occluded  while  its  uncertainty  is 
moderate.  EKFSLAM  and  FastSLAM  have  almost  similar  ranging  accuracy;  however,  EKF  have 
slightly  better  performance  in  NLOS  areas.  EKFSLAM  allows  FP  locations’  uncertainty  to  be 
remembered  as  it  saves  and  keeps  the  track  of  their  information  in  terms  of  the  covariance 
matrix.  Moreover,  the  correlation  between  states  in  the  covariance  matrix  allows  the  FP  location 
estimate  to  be  updated  in  its  NLOS  area.  Whenever  they  reappear,  this  information  is  re-called 
and  the  estimate  can  also  be  highly  improved  due  to  its  correlation  with  MN  location  and  other 
FP  estimates  specially  if  other  FP  locations’  estimates  have  obtained  good  accuracy.  However, 
in  the  FastSLAM  method,  whenever  a  FP  goes  out  of  sight,  the  resampling  may  discard  most  of 
the  particles  relating  to  that  FP,  so  eventually  the  entire  FP  location  history  map  is  lost  and 
cannot  be  recovered  if  the  FP  reappears. 
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Figure  0-31  MN-ANs  actual  and  estimated  distances  for  case  1 


The  position  error  bound  (PEB),  obtained  by  trace  of  the  corresponding  submatrix  of  the 
BFIM  can  be  utilized  to  evaluate  the  localization  result  error  instead  of  the  whole  BFIM  which 
also  includes  nuisance  parameters  information.  As  an  extension  of  Eq.  (0-113)  for  OWLS,  we 
modified  this  bound  for  both  the  MN  location  estimate  and  FP  locations  estimates,  given  by 


:{(x,-x,)2+(^-v,)2}>rrj[/(qi:,)^ 


MN 

2x2 


0-139 


E  {m,,x  -  )"  +  (mi,y  -  ™i,y,t  M  -  H  J  H:, 


FP,i 

2x2 


where 


-| MN  r  ■ 

ar|d  Jtot  (q1:r) 

J2x2  L 


FPj 


are  the  2x2  submatrixes  of  inverse  FIM 


2x2 


corresponding  to  the  MN  location  and  / th  FP  location  variables,  respectively.  Since  the  true  joint 
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probability  density  function  is  not  available;  to  calculate  the  left  side,  the  expectation  is 
approximated  by  averaging  over  N=50  Monte  Carlo  simulation  runs, 


{(**  -  **  )2 + (y*  ~  y<  )2 } *  -L  E  ((**  -  **  )2 + (yt  ~  P,  f ) 

’  n= I  n 


0-140 


The  50-run  average  position  error  plots  for  the  MN  location  and  FP  locations  are  shown 
in  Figure  0-30  for  both  EKF-SLAM  and  FastSLAM  estimates.  It  can  be  seen  that  FastSLAM, 
unlike  its  overconfident  uncertainty,  outperforms  EKFSLAM  in  terms  of  MN  location  error, 
except  in  interval  [80,100]  when  the  information  from  re-seen  AN01  and  AN02  must  be  restored. 
In  loop  closure,  when  the  MN  moves  through  unknown  trajectory,  and  at  some  points 
encounters  FPs  which  are  not  seen  for  a  long  time,  maintaining  the  correlation  between 
estimates  is  very  essential  in  a  SLAM  algorithm.  In  EKFSLAM,  it  is  maintained  directly  in  the 
covariance  matrix  while  in  FastSLAM,  it  is  maintained  through  particles  diversity.  Unfortunately, 
in  FastSLAM,  new  observations  cannot  improve  the  FP  location  estimates  observed  prior  to  this 
point,  thus  their  correlation  to  the  particles  degrades  and  corresponding  particles  are  discarded 
eventually  through  resampling.  Moreover,  PEB  plots  can  reveal  how  much  new  information  can 
improve  the  localization  accuracy  and  if  it  worth  to  use  it.  Comparing  PEB  results  for  4-FP  and  3 
FP-  1  AP  scenarios  in  Figure  0-32,  it  can  be  inferred  that  for  the  given  example,  the  knowledge 
of  a  single  FP  location  improves  position  error  MN  location  estimate  approx,  by  0.2  meters. 


Time  Step 


a)  4  FP 


Time  Step 

b)  3  FPs-1  AP 


Figure  0-32  MN  position  error  in  meters  for  case  1 
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Figure  0-33  shows  the  FP04  position  error  and  PEB  obtained  from  the  BFIM.  It  is 
illustrated  that  in  a  Gaussian  linear  system,  FP  PEB  will  converge  to  zero;  however,  as 
FastSLAM  or  EKFSLAM  cannot  guarantee  this  convergence  due  to  system  nonlinearity. 

Another  interesting  observation  from  comparison  between  Figure  0-33  and  Figure  0-32,  both 
EKFSLAM  and  FastSLAM  can  achieve  PEB  for  the  MN  location  with  good  approximation  but 
not  for  the  FP  location  estimates.  Therefore,  the  total  Euclidian  distance  between  true  state 
vector  and  estimated  state  vector  is  dominated  by  the  FP  location  error  as  shown  in  Figure 
0-31.  It  must  be  mentioned  that  effect  of  initial  conditions  is  not  considered  in  MC  simulation  and 
the  results  is  averaged  for  50  realizations  of  measurement  model  with  constant  initial  condition 
for  MN/FP  locations. 


Time  Step 


Time  Step 


a)  4  FPs  b)  3  FPs-  1  AP 

Figure  0-33  FP04  position  error  in  meters  for  case  1 
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Time  Step 


a)  4  FPs 


Time  Step 

b)  3  FPs-  1AP 


Figure  0-34  Total  MN/FPs  location  error  in  meters  for  case  1 


Another  measure  to  characterize  the  SLAM  performances  is  known  by  NEES  s  ,  earlier 
outlined  in  Eq.  (0-114).  Nonlinearity  nature  of  SLAM  prevents  it  to  find  an  efficient  solution  by 
FastSLAM  or  EKFSLAM.  However,  averaged  NEES,  s  ,  as  defined  in  Eq.  (0-115)  over  N  runs 

can  still  be  used  to  evaluate  the  SLAM  performance  for  inconsistency.  It  is  due  to  the  fact  that 
since  the  error  sequence  is  correlated,  a  single  run  NEES  cannot  follow  the  chi-square 
distribution.  Given  N=50  runs,  the  consistency  of  both  solutions  are  evaluated  by  averaged 
NEES  as  given  in  Eq.  (0-114)  and  (0-115)  for  the  MN  location,  ANs  location,  and  the  whole 
state  vector.  Given  the  hypothesis  of  a  consistent  linear-Gaussian  filter,  averaged  NEES  has  a 

chi-squared  density  with  ,/Vdim(q,)  degrees  of  freedom 

Thus,  with  N  =50,  the  two  sided  95%  probability  confidence  region  for  the  2-dimensional 
location  (for  the  MN  or  any  FP),  is  bounded  by  the  interval  [1 .48,2.59],  for  the  2  +  2  N^- 

dimensional  state  vector  is  bound  by  the  interval  [6.93,9.15].  If  e  rises  significantly  higher  than 

the  upper  bound,  the  filter  is  optimistic,  if  it  tends  below  the  lower  bound,  the  filter  is 
conservative.  As  shown  in  Figure  0-35  and  Figure  0-36,  there  is  a  discernible  difference 
between  EKFSLAM  and  FastSLAM  in  terms  of  NEES.  The  measure  of  performance  obtained  by 
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EKFSLAM  for  MN  location  and  FP  locations,  corresponding  marginal  posterior  PDF  covariance 
matrix,  is  far  superior  since  the  calculated  NEES  stays  within  consistency  bounds  and  close  to 
the  dimensionality  of  state  space  (2,  here  for  the  FP  location  or  MN  location).  NEES  statistics 
indicated  that  FastSLAM  underestimates  the  covariance  matrix  due  to  loss  of  particles  diversity. 
However,  EKFSLAM  performance  is  consistent  in  state  space  of  MN  location  and  FP  location 
separately;  it  underestimates  the  posterior  covariance  matrix  for  overall  state  space  as  shown  in 
Figure  0-37. 


Time  Step 


Figure  0-35  MN  location  averaged  NEES  from  EKFSLAM  and  FastSLAM  for  case  1 
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Figure  0-36  FP04  location  NEES  from  EKFSLAM  and  FastSLAM  for  case  1 
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Figure  0-37  Total  state  NEES  from  EKFSLAM  and  FastSLAM  for  case  1 


System  model  2:  synchronized  ANs,  MN  not  synchronized  with  ANs,  known  sight 
conditions 

The  previous  case  modeled  an  idealistic  SLAM  problem  where  the  MN-ANs  links  are 
synchronized.  In  practical  situation,  the  MN  clock  is  subject  to  drift,  leading  to  bias  error  in  TOA 
or  range  measurements.  As  earlier  explained  in  chapter  03,  the  random  effect  of  MN  clock  drift 
is  modeled  as  a  random  walk  process  which  can  be  estimated  to  mitigate  its  effect  on  the  MN 
location  estimate  and  FP  mapping.  The  effect  of  new  augmented  nuisance  variable  is 
investigated  and  the  result  is  compared  with  the  case  where  it  is  ignored. 

EKFSLAM  solution  for  case  2 

For  the  expanded  state  vector,  q  =  [p,  b,  m/;Jr ,  the  Gaussian  posterior  PDF  is 


approximated  by 
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where  p(|M  and  !Lpp  are  the  mean  and  covariance  matrix  of  apriori  PDF,  stated  by 
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Dt=HTllt_,HT+R, 

K,  =  HtD~x 

H  is  the  Jacobian  matrix  of  range  measurement  function,  evaluated  at 
q  =[p,M  b,M  m/>(_1]r ,  where  the  elements  are  defined  in  Eq.  (0-100)  and  (0-101). 

FastSLAM  solution  for  case  2 

Since  the  range  offset  is  common  between  all  ANs  measurement,  the  joint  posterior  PDF 
is  factorized  based  on  Rao-Blackwellizied  theorem  for  case  2  as  follows 

NFP  0-143 

P(p,t,blt,m/p  |  Zj  =  p(pit,blt  |  ZVt)Y[p(mfP,,  I  P.:rAr>Zl:,) 

!=1 

The  trajectory  of  MN  locations  and  range  offsets  is  represented  by  weightening  samples 
and  the  map  is  computed  analytically  using  the  EKF  algorithm.  The  mth  particle  of  the  posterior 
PDF,  at  time  step  t ,  represents  not  only  the  history  of  MN  locations  but  also  the  history  on 
range  offsets. 

The  apriori  PDF  is  approximated  by  NP  independent  and  identically  distributed  (i.i.d) 


samples  (particles) 


p(tux  i  M4,. -*"'*) 

7V  o  m= 1 


0-144 
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pin  b  |p[m]  b[m]  m[",]) 

which  is  drawn  from  proposal  distribution,  ' irM’  '  1  ’  &  .  Then  similar  to  previous 

case,  for  each  particle,  an  EKF  update  is  performed  over  the  observed  ANs  as  a  simple 
mapping  operation  under  the  known  MN  location  and  range  offset  assumption. 

For  the  same  trajectory  used  for  case  1,  the  proposed  EKFSLAM  and  FastSLAM 
solution  are  applied  for  a  single  Monte  Carlo  simulation  run  of  system  model  1 .  The  initial 
conditions,  measurement  noise,  and  process  noise  parameters  have  kept  the  same  as  the  ones 
used  in  case  01  for  comparison  purposes,  except  the  measurement  with  range  offset  model 
according  to  Eq.  (0-67)  is  used  for  generating  range  measurement  data. 

Figure  0-38  (a)-(d)  compare  the  EKFSLAM  and  FastSLAM  performance  in  a  single  run 
MC  for  4  FPs  and  3  FPs-1  AP  scenarios  where  the  range  offset  is  tracked  along  with  other 
location  parameters,  and  Figure  0-38  (e)-(f)  illustrates  their  performance  where  the  range  offset 
effect  is  not  ignored.  For  the  4  FPs  scenario,  the  FastSLAM  shows  a  slightly  better  performance 
than  the  EKFSLAM,  as  shown  in  Figure  0-38  (a)-(b),  while  their  performance  is  comparable  in 
the  3  FPs- 1  AP  scenario,  as  shown  in  Figure  0-38(c)-(d).  Since  the  FastSLAM  tracks  its 
solutions  in  a  wide  area  of  room  through  particles,  it  shows  higher  robustness  to  ambiguous  and 
uncertain  situation  in  comparison  to  the  EKFSLAM  such  as  4  FPs  case  when  there  is  no 
knowledge  about  source  locations.  The  effect  of  range  offset  is  observed  in  Figure  0-38(e)-(f), 
where  MN  localization  error  is  increased  to  more  than  1  meter  for  most  of  the  MN  trajectory 
while  it  is  less  than  1  meter  in  average  when  the  range  offset  is  modeled,  as  illustrated  Figure 
0-38(a)-(b).  The  range  offset  estimates  along  trajectory  for  both  EKFSLAM  and  FastSLAM  are 
shown  in  Figure  0-37. 
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EKFSLAM  with  0  Anchor  Points, 4  Feature  Points 


FastSLAM  with  0  Anchor  Points, 4  Feature  Points 
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Figure  0-38  Example  of  MN  tracking  and  ANs  mapping  for  system  model  2 
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a) 


Time  Step 

c) 


Time  Step 


b) 


Time  Step 

d) 


Figure  0-39  Bias  tracking  for  case  2 


For  case  model  2,  Figure  0-40  shows  the  actual  MN-ANs  distances  and  the 
corresponding  estimates  obtained  by  EKFSLAM  and  FastSLAM.  For  the  NLOS  area,  the 
ranging  information  is  blocked  with  a  grey  rectangular.  It  is  shown  range  estimation  is  exact 
(with  error  less  than  0.5  meter)  for  all  methods  even  when  bias  effect  is  ignored,  while  the  MN 
localization  error  reaches  to  2  meters. 
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Figure  0-40  MN-ANs  actual  and  estimated  distances  for  case  2 


The  50-run  average  position  errors  for  the  MN  location,  FP  04  location  and  total  state  is 
shown  in  Figure  0-41,  Figure  0-42,  Figure  0-43,  for  both  EKF-SLAM  and  FastSLAM  estimates. 
Unlike  the  expectation  from  single  MC  run  results,  MN  position  error  obtained  by  FastSLAM  is 
higher  than  position  error  obtained  by  EKFSLAM.  This  is  because  FastSLAM  performance  is  not 
consistent,  as  also  observed  in  case  1,  due  to  particle  diversity  loss.  FastSLAM  particles  does 
not  represent  a  momentary  MN  location  and/or  range  offset;  they  represent  the  entire  MN 
trajectory  and  range  offset  history.  The  dimension  of  state  space  increases  with  time  while  the 
number  of  particles  representing  this  large  state-space  dimension  become  insufficient.  This  is 
observed  as  the  degradation  of  FastSLAM  performance  in  the  case  2  in  comparison  with  its 
performance  in  the  case  1  where  it  outperformed  EKFSLAM  in  term  of  MN  location  error 
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(excluding  loop  closure  areas).  The  results  comply  with  consistency  test  obtained  by  NEES  test 
shown  in  Figure  0-44,  Figure  0-45,  and  Figure  0-46. 
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Figure  0-41  MN  position  error  in  meters  for  case  2 


Time  Step 

Figure  0-42  FP04  position  error  in  meters  for  case  2 
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Figure  0-43  Total  state  error  in  meters  for  case  2 
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Figure  0-44  MN  location  averaged  NEES  from  EKFSLAM  and  FastSLAM  for  case  02 


121/154 


10 


Single  FP  NEES  (FP04),  4FP  case  02 


95°/ 

upper  bound 

lower  bound 

SLAM,  4FP 
TSLAM,  4FP 

- EKF 

- FAS 

— - — 

0  20  40  60  80  100 

Time  Step 

Figure  0-45  FP04  location  NEES  from  EKFSLAM  and  FastSLAM  for  case  02 
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Figure  0-46  Total  state  NEES  from  EKFSLAM  and  FastSLAM 


System  model  3:  unsynchronized  ANs,  MN  not  synchronized  with  ANs,  known  sight 
condition 

In  the  third  case,  the  system  model  is  extended  to  when  the  MN  free  running  clock  is  not 
synchronized  with  ANs,  and  ANs  are  not  time  synchronized  with  each  other.  To  model  the  range 
measurement  data,  a  single  bias  for  each  MN-AN  link  can  model  both  MN  and  ANs  clock  drifts. 

As  earlier  described  in  Section  0,  the  MN  location  state  p,  is  modeled  as  a  first-order  Markov 
process  ruled  by  P,=Pm+v„  where  the  \t  Q  N(0,Q  )  .  Based  on  smooth  trajectory 
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assumption,  the  range  offset  vector,  b,  =[V  bN^  t  J  can  also  be  modeled  as  Nan 
independent  first-order  Markov  processes,  given  by  b,  =  b,_,  +vA(,  where  ybt  □  N(0,Qb )  and 
Qb  =  Jwg  {ofA  •  •  ■  cr^v  A  | .  The  range  offset  evolves  independently  from  MN  location  states. 
The  current  state  vector  of  unknown  variables  is  denoted  by  q,  =  [p,  b;  m/p]  ,  where  m  //; 


is  the  stacked  vector  of  the  FP  locations,  defined  by  [in.,  m,r  ■■■  m,,  T. 

3  L  nap+ 1’x  nap* l>y  NAN’X  Nan,v} 

Modified  version  of  EKF-SLAM  and  FastSLAM  as  suboptimal  solutions  for  this  system  model 
are  proposed  based  on  reasonable  assumptions  that  allow  the  evaluation  of  posterior  PDF  with 
efficient  computation. 

EKF-SLAM  solution  for  case  3 

For  the  defined  state  vector,  q  =[p,  b,  m/;Jr ,  the  EKFSLAM  approximates  the 


posterior  PDF  by  multivariate  Gaussian  PDF  where  it  means  and  covariance  are  given  by 
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0-147 


F  = 


10  0  0 
0  0  0 
0  0  10 


2+Nan  +2Nfp 


Q  =  diag  \Qr  aj 
v,  =  hi.„_,ht+r, 

K,  =  htd;‘ 

H  is  the  Jacobian  of  range  measurement  function,  evaluated  at  q  =  [p,M  b 


t\t-l 


m^-i] 


previously  defined  in  Eq.  (0-103)  and  (0-104). 

FastSLAM  solution  for  case  3 

The  traditional  FastSLAM  method  is  based  RB  factorization  where  the  EKFs  used  for 
FPs  location  is  only  conditioned  on  MN  trajectory  which  is  represented  by  particles.  The 
simulation  results  from  case  02  showed  that  FastSLAM  performance  is  very  sensitive  to  the 
dimension  of  state  space  which  is  represented  by  particles.  The  FastSLAM  degradation  due  to 
loss  of  particle  diversity  is  observed  for  case  3  where  the  particle  space’s  dimension  is 

increased  to  Nan  xt  +2 xt  where  is  the  dimension  of  range  offset  sequence  b1:,  and  2xt  is 


the  dimension  MN  trajectory  p1:, .  In  this  regard,  we  have  modified  the  assumption  of  FastSLAM 

in  such  a  way  that  the  particle  dimension  can  stays  the  same  as  original  version  (as  in  case  1), 
2  xt ,  but  the  conditioned  EKF  part  of  FastSLAM  is  modified  to  jointly  estimates  the  low 
dimension  FP  location  and  its  range  offset  recursively.  Therefore,  The  Rao-Blackwellization 
factorization  for  joint  posterior  PDF  is  modified  as 

Nan  0-148 

p(p,„b,nmJp\zll)  =  piptjzj  p(mfp.,biht\p,t,zht) 

i=NAP+ 1 

The  m th  particle  of  posterior  PDF,  at  time  t ,  represents  the  history  of  MN  locations 
denoted  as  pLj"] ,  which  additionally  contains  the  Gaussian  parameters  (the  mean  and 
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covariance  matrix)  of  Ay’1 ,2[”]  ,  [  of  jointly  FPs  location  and  its  range  offset 


yFP’1  iy  FP 


w*»Pw 

First,  the  FastSLAM  extend  the  MN  trajectory  by  drawing  new  sample  according  to  the 
motion  statistics  as 

r  1  r  .  0-149 

p,1"1  □  Pi  p,  i  p,-,1"1) 

Next,  for  each  particle,  an  EKF  is  implemented  to  update  FP  estimates  and  range 
offsets.  The  output  of  this  stage  is  ,  defined  by 


yl>]  _  yl>]  _|_  J£l">]£)[m]j£[m],T 

where,  for  observed  FPs,  i  =  Nap+1,...,Nan , 


0-150 


U[m]  =  Ll[m] 


0-151 


VLmJ  _  VL»'J  +  tT 2  p2 
+Cri,4-e3x3 


pM  _  y[»'l 

for  observed  APs,  /  =  l,...,iV 


AP 


m  m 

Ml,  =  M\ A 

yM  _  y['«]  i  _2 

^ i,t  ~  ^i.t-1  +  ai,b 

and 

P)lm]  _  zj['»]  y[m] 

11 1  t\t-\ 11t  ‘  ^i,z 
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Next,  the  samples  weights  w|"'J  are  calculated  according  to  importance  function  which  is 
based  on  the  conditional  probability  of  p(zt  |p[m]).  Finally,  resampling  is  accomplished  with 

replacement  from  the  set  jpo™1}^ ,  including  their  associated  AN  locations’  and  range  offsets’ 


estimates  considering  each  particle  has  probability  proportional  to  w\ 
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For  the  same  trajectory  used  for  case  1  and  2,  the  proposed  EKFSLAM  and  FastSLAM 
solutions  are  applied  for  a  single  Monte  Carlo  simulation  run  of  system  model  03.  The  initial 
conditions,  measurement  noise,  and  process  noise  are  similar  to  case  01  for  comparison 
purposes,  except  the  measurement  data  is  generated  based  on  the  system  model  3. 

Figure  0-45  (a)-(d)  compare  the  EKFSLAM  and  FastSLAM  performance  for  a  single  run 
MC  for  4  FPs,  and  3  FPs-1  AP  scenarios  where  the  range  offsets  are  tracked  along  with  other 
location  parameters;  whereas  Figure  0-47  (e)-(f)  illustrates  their  performance  where  the  range 
offset  effect  is  ignored.  For  4  FPs  scenario,  modified  FastSLAM  and  EKFSLAM  shows  a  similar 
performance,  as  shown  in  Figure  0-47  (a)-(b),  while  the  EKFSLAM  have  slightly  better 
performance  for  3  FPs-  1AP  ,  as  shown  in  Figure  0-47  (c)-(d).  The  effect  of  range  offset  is 
observed  in  Figure  0-47  (e)-(f);  MN  localization  error  increases  to  4  meters,  EKFSLAM  and 
modified  FastSLAM  which  model  and  track  range  offsets  provide  less  than  1  meter  localization 
error.  The  range  offset  estimates  within  trajectory  for  both  EKFSLAM  and  FastSLAM  is  shown  in 
Figure  0-48. 
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Figure  0-47  Example  of  MN  tracking  and  ANs  mapping  for  system  model  3 
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c)  EKFSFAM,  3  FPs-  1AP 


FASTSLAM,  4FP  case  03 


b)  Modified  FastSFAM,  4  FPs 
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d)  FastSFAM,  3  FPs-  1AP 


Figure  0-48  Bias  tracking  for  case  3 


For  case  3,  Figure  0-49  shows  the  actual  MN-ANs  distances  and  the  corresponding 
estimates  obtained  by  EKFSLAM  and  FastSLAM.  It  is  shown  range  estimation  is  exact  (with 
error  less  than  0.5  meter)  for  all  methods  even  when  range  offset  is  ignored.  Both  EKFSLAM 
and  FastSLAM  have  comparable  performance  in  ranging  accuracy. 

The  50-run  average  position  error  results  for  the  MN  location,  FP  04  location  and  total 
state  vectors  are  shown  in  Figure  0-50,  Figure  0-51,  and  Figure  0-52,  for  both  EKF-SLAM  and 
FastSLAM  estimates.  MN  position  errors  obtained  by  both  the  FastSLAM  and  EKFSLAM  are 
comparable  and  approximately  near  the  optimal  PEB  bound  obtained  by  BFIM,  except  loop 
closure  area  where  the  FastSLAM  fails.  In  terms  of  FP  position  error,  the  EKFSLAM  not  only 
outperforms  FastSLAM  but  also  can  reach  the  PEB.  NEES  results  obtained  for  the  modified 
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FastSLAM  and  EKFSLAM  shown  in  Figure  0-53, Figure  0-54, Figure  0-55  evaluates  the 
consistency  of  both  methods. 
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Figure  0-49  MN-ANs  actual  and  estimated  distances  for  case  3 


The  FastSLAM  performance  is  not  consistent  because  of  the  same  reasons  already 
explained  for  case  1  and  2.  On  the  other  hand,  EKFSLAM  consistency  is  degraded  due  to 
expansion  of  state  space  in  comparison  with  previous  cases,  where  the  multivariate  Gaussian 
PDF  will  no  longer  can  represent  joint  MN  location,  FP  location  and  AN  range  offsets  state 
space.  The  uncertainty  obtained  from  EKFSLAM  is  underestimated  since  the  NEES  is  near  the 
lower  bound  for  FP  04  location  and  completely  under  lower  bound  for  overall  state  space. 
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Time  Step 

Figure  0-50  MN  position  error  in  meters  for  case  3 


Time  Step 

Figure  0-51  FP04  position  error  in  meters  for  case  3 
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Figure  0-52  Total  state  error  in  meters  for  case  3 
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Figure  0-53  MN  location  NEES  from  EKFSLAM  and  modified  FastSLAM  for  case  3 
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Figure  0-54  FP04  NEES  from  EKFSLAM  and  modified  FastSLAM  for  case  3 
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Figure  0-55  Total  state  NEES  from  EKFSLAM  and  modified  FastSLAM  for  case  3 
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System  model  4:  unsynchronized  ANs,  MN  not  synchronized  with  ANs,  unknown  sight 
conditions 

In  the  system  model  4,  it  is  assumed  that  the  MN  receives  signal  of  opportunity  from 
unsynchronized  network  in  a  non-stationary  mixed  LOS/NLOS  propagation.  The  statistic  of  the 

observables  zt  changes  over  the  time  due  to  random  variation  of  clock  drift  and  the  AN  sight 

conditions.  In  such  a  scenario,  optimal  tracking  algorithm  requires  the  definition  of  additional 
dynamic  models  that  describe  the  evolution  of  range  bias  as  well  as  ANs  sight  state.  This  leads 
to  a  composite  dynamic  model  where  the  state  is  composed  of  the  joint  MN  position,  AN 

location  ,  range  bias  and  sight  state  variables  q,  =  [p,  b,  s,  m]r  where 

s(=[su  s2t  ...  V,,v*]  is  the  vector  of  ANs’  sight  states  at  time  step  t. 

This  composite  localization  system  is  described  by  a  nonlinear  jump  Markov  system 
where  the  state  and/or  the  measurement  model  depend  on  a  driving  Markov  chain.  Following 
our  assumptions,  in  the  particular  case  herein  considered,  only  the  range  measurement  model 

z,  depends  on  the  jumping  feature  s, .  In  fact,  the  PDF  of  measurement  z,  is  driven  by  a  the 
discrete  process  s, ,  while  the  state  p,  is  assumed  to  be  independent  of  s, .  Once  the  JMS  is 
defined,  the  evolution  of  joint  state  q,  can  be  tracked  by  a  Bayesian  filter. 

As  outlined  in  Section  0  for  the  system  model  4,  the  MN  location  p,  is  modeled  as  a 
first-order  Markov  process  ruled  by  Pf  “Pm  +v<>  where  v,  □  N(0,Q  ) .  Based  on  the  smooth 

trajectory  assumption,  the  range  offset  vector,  b,  =[/>,,  ^VIA,«]  is  described  by  Nan 

independent  first-order  Markov  processes,  given  by  b,  =  b,_,  +  vJ>f ,  where  \b  t  □  N(0,Qb)  and 
Qb  =  diag  jcr,2/;  •  •  •  <j2n ^  b  J .  The  range  offset  evolves  independently  from  MN  location  states. 
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In  the  same  manner,  each  AN  sight  state  variable  s;,  is  modeled  as  a  binary  Markov  chain.  The 
Markov  chain  process  is  described  by  a  2x2  transition  probability  matrix  n  which  is 
completely  defined  by  the  transition  probabilities  P(si  t  =  0 1  sit_x  =  0 )-  p0  and 

P(sj  t  =  1 1  sf  M  =  1)  =  px .  Under  the  independence  assumption  between  the  MN-AN  links,  the 
whole  sight  state  vector  s,  is  determined  by  a  first-order  Markov  chain  with  transition 

Nan 

probabilities  P( s,  =  c  |  sM  =  d)  =  =  c/ 1  -  Vi  =  di)  for  c  =  Lci  cnan \  > 

i= 1 

d  =  ^dx  dNm  J  where  di , ci  c  {0, 1} . 

The  overall  state  vector  of  unknown  variables  is  denoted  by  q,  =  [p,  b,  s,  m^]  , 
where  m  fp  is  the  stacked  vector  of  the  FP  locations,  defined  by 

m,,  ■■■  m.r  1  .  The  transition  PDF  from  previous  state 

L  iy AP+  l’x  iyAP+\->y  iyAN'X  iyAN*y 

q,_i=[p,_i  b,_,  s,_,  m]r  to  the  next  one  q,  =  [p,  b,  s,  m]r  is  given  by 

p{(\t  |  qM)  =  p( p,  |  p,_|)/?(br  |  b,_x) /?(s,  |  sM)  under  the  assumption  that  MN  trajectory  evolves 
independently  from  both  clock  drift  and  ANs’  sight  state  processes.  The  JMS  state  q,  is  hidden 
into  the  A^-link  observation  vector  z,  composed  by  the  conditionally  independent 

Nan 

measurements  p( z,  |  q()  =  ]^[ p(zi  t  |  q,)  where  each  p(zit  |  q,)  has  to  be  evaluated  according 

i=i 

to  the  observation  model  defined  in  Eq.  (0-70)  and  (0-71) 

An  optimal  estimate  of  q,  can  be  extracted  from  the  whole  set  of  available  measurement 

by  timer,  zl7  □  {z,  z2  ...  zj  where  zA  ={  zlk  z2k  ...  zNjl/k}  denotes  range 

measurements  at  time  step  k  from  all  ANs  in  NLOS  or  LOS  with  the  MN,  where  the  AN  sight 
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conditions  are  unknown  to  the  MN.  The  estimate  is  obtained  by  evaluating  the  a  posterior  pdf 
p( qf  |  zH)  which  can  be  expressed  as  a 


P(Pr™fp  I  zi,)  =  ‘’mfP  I  SL>Z1  ;)PK  I  Zl:r) 

1=1 

where  p(s[.t  |  zl7)  is  the  probability  of  a  particular  sight  state  sequence,  sjt ,  given  the 
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measurement.  Ns  =  2Nan  is  the  total  number  of  all  possible  ANs  sight  state  combinations  at  each 
time  step.  Each  possible  ANs  sight  state  combination  is  also  referred  as  sight  mode.  As  earlier 
discussed,  the  optimal  solution  for  p(pt,vafp  |  zl7)  cannot  be  found  analytically  in  a  closed  form 

due  to  nonlinear  feature  of  measurements.  Moreover,  the  number  of  mixture  components  in  the 
PDF  sum  of  Eq.  (0-153)  grows  exponentially  with  time.  Hence,  a  practical  solution  has  to  be 
suboptimal  Following,  two  solutions  for  this  jump  Markov  SLAM  problem  based  on  EKF  and  PF 
are  proposed. 

IMM-EKFSLAM 

A  practical  implementation  of  JMS  solution  is  to  limit  the  growth  of  number  of  mixture 
components  by  merging  of  the  sight  state  mixture  components.  This  type  of  algorithm  is 
generally  known  as  interactive  multiple  model  (IMM)  estimator.  The  proposed  IMM-EKFSLAM 
algorithm  employs  the  EKF  to  linearize  the  nonlinear  system  and  filter  every  mode  and  then 
uses  the  nonlinear  IMM  algorithm  to  fuse  the  estimated  state  of  the  filtering  for  each  mode,  and 
finally  to  get  an  overall  estimated  state.  We  define  the  state  vector  for  each  sight  mode  as 

qs  /  =  [p,  b,  m]r  ■  Before  filtering  for  each  sight  mode,  the  IMM  obtains  a  hybrid  prior 

estimate  by  mixing  of  modes  where  the  mean  and  covariance  of  Gaussian  posterior  PDF  at  time 
step  t-\  for  sight  mode  j  are  defined  by 
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Ns 


aw  =  y  u^j)  a(i) 

i=l 

Ns 

yO)  _y  r y(0  ,(£(0  ~au)  ¥a(0  -nU)  ')rl 

;=i 

It  is  easier  to  address  sight  modes  by  an  integer  number  instead  of  a  binary  vector.  In  this 
regard,  we  define  rt  a  {\,...,NS }  as  an  integer  index  for  the  binary  vector  of  ANs’  sight  state  s, . 

With  this  definition,  the  weights  as  the  mixing  probabilities  are  given  by 


p(r,_t=i\ri_,=  j,ztl)  = 


n  ,x; 
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2XX:; 

i= 1 

where  n  =  n1  ®...0nJV"r  and  n\]_\  is  the  y-th  mode  probability  which  obtained  after  filtering  at 

previous  time  step.  The  next  phase  is  extended  Kalman  filtering  based  on  each  sight  mode, 
where  the  prediction  and  update  stages  are  stated  as 
Filter  predict  stage: 


q«  = 
4//-1 


’  P,M  ’ 

(0 

P/-UT-1 

Vi 

= 

3- 

T 

i _ 
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3» 

T 

i _ 

s<0,M=s!'V,+^re^ 

Filter  update  stage: 


D f=HU)  I(V1(^(0)r+^(° 


0-157 


K 


(0 


’  P*  " 

(0 

i 

T 

_ i 

a(i) 

4S,^= 

= 

1 

<E 

_ 1 

i 

3> 

T 

i _ 

-|(0 


+  K"'[(Z,-A(Cm)] 


Finally,  the  output  of  IMM-EKFSLAM  filter  is  calculated  as  the  mean  and  the  covariance 


of  the  Gaussian  mixture  of  a  mode-based  posterior  PDF  as  follows 
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Ns 


q,,k  = 


i= 1 
N~ 


2*  =  Z^(i)  P'V  +(qU  -i^XqU  -<w)r] 

/=1 

where  the  mode  probabilities  are  calculate  using  the  Bayes  rule  as 


N . 


«u>  □ 


A/Zn, 


0-759 


7=1 


2A/Zn.v«"! 


7=1  1=1 


where  A)  =  p(zt  |  z Vt,rt  =  /)  □  JV((zf  -h{ q*'jM),D^  ) .  Moreover,  the  sight  mode  is  estimated  by 


(O'! 


rt  =  arg  max 
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IMM-FASTSLAM 


The  IMM-FastSLAM  is  proposed  as  an  extension  of  jump  Markov  PF  (also  known  as 
MM-PF)  for  SLAM  problem  with  switching  dynamic  models.  However,  the  RB  factorization  for 
posterior  PDF  is  modified  by 


' AN 

I  2l:t  )  —  P(Po.t  ’  *0:f  I  Z1 :()  h[  P  (®  fp  ,i  ’  ^ '/  ,0:t  I  Po:?  ’  ^0:t  ’  Zl:l  ) 

i=^AP+ 1 

In  first  stage,  the  algorithm  generates  a  random  set  {s[m]}AV  based  on 

W  V  J  m= 1  v  1  '  777=1 


0-161 

and  the 


state  mode  transitional  probability  matrix  n .  Next,  IMM-Fast  SLAM  performs  a  PF  filtering 
conditioned  by  mode  particles.  For  each  particle,  the  EKF  is  applied  to  update  the  estimation  of 
range  offsets  and  FP  locations.  Finally,  particle  weights  are  updated  based  on  EKF  estimates 
and  the  resampling  is  done  as  before.  The  state  vector  can  be  calculated  through  the  weighted 
means: 


136/154 


0-162 


m=l 


m= 1 


m= 1 


where  the  FPs  location  estimates  and  range  offset  estimate  bjt[m]  is  obtained  from  the 

mean  of  the  EKF  solution  for  m th  particle. 

A  similar  MN  trajectory  and  ANs  map  used  for  the  previous  case  is  considered  for 
evaluating  the  proposed  methods  for  case  4.  The  initial  conditions,  measurements  and  transition 
model  parameters  are  kept  the  same  as  before  except  the  measurement  data  is  generated 
according  to  Eg.  (0-70)  and  (0-71)  where 
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Figure  0-56(a)-(b)  compare  the  EKFSLAM  and  modified  FastSLAM  performance  for  a 
single  run  MC  for  4  FPs  scenario  when  the  ANs  sight  mode  is  detected  by  local  maximum 
likelihood  (LML)  criteria 

r,=argmaxA'=argmax/>(z?|r?  =/,z1:M)  °'16 


where  z'  is  range  estimate  obtained  based  on  i  -th  sight  mode.  Then  MN  location  and  FP 

locations  estimates  are  chosen  based  on  the  detected  sight  mode.  For  the  same  measurement 
data  set,  the  proposed  IMM-EKFSLAM  and  IMM-FastSLAM  are  implemented.  It  can  be 
observed  that  localization  accuracy  is  improved  for  both  in  comparison  with  LML-EKFSLAM  and 
LML-FastSLAM.  In  particular,  the  confidence  region  obtained  by  IMM-EKFSLAM  reaches  the 
confidence  region  obtained  by  BFIM. 
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EKFSLAM  with  0  Anchor  Points, 4  Feature  Points 
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IMM-EKFSLAM  with  0  Anchor  Points, 4  Feature  Points 
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FastSLAM  with  0  Anchor  Points, 4  Feature  Points 
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IMM-FastSLAM  with  0  Anchor  Points, 4  Feature  Points 
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Figure  0-56  Example  of  MN  tracking  and  ANs  mapping  for  system  model  4 
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Figure  0-54  (e)-(f)  illustrates  the  performance  of  EKFSLAM  and  FastSLAM  (as  explained 
in  case  1)  where  the  range  offset  effect  is  not  ignored  and  the  MN  has  no  knowledge  about  ANs 
sight  conditions.  It  can  be  seen  that  for  the  area  that  ANs  are  shadowed  by  the  wall,  bias  delay 
effect  leads  to  large  MN  localization  error.  For  this  MC  simulation,  PF-based  methods  (IMM- 
FastSLAM  and  LML-FastSLAM)  have  better  performance  in  terms  of  MN  localization  accuracy; 
however,  EKFSLAM  outperforms  in  terms  of  FP  locations  and  their  uncertainty  estimation. 

Actual  ANs  sight  state  and  their  estimates  from  LML-EKFSLAM,  LML-FastSLAM,  IMM- 
EKFSLAM,  and  IMM-FastSLAM  are  depicted  in  Figure  0-57 ,  Figure  0-58  .Figure  0-59 ,  and  Figure 
0-60.  It  is  shown  that  IMM-based  methods  outperform  LML-based  methods  since  they  can  take 
advantage  of  Markovian  process  between  sight  states  to  incorporate  the  information  from 
previous  measurements. 
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Figure  0-57  Sight  state  tracking  according  to  LML-EKFSLAM 
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Figure  0-58  Sight  state  tracking  according  to  LML-FastSLAM 
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Figure  0-59  Sight  state  tracking  according  to  IMM-EKFSLAM 
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Figure  0-60  Sight  state  tracking  according  to  IMM-FastSLAM 
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Figure  0-61  MN  position  error  in  meters  for  case  4 
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Figure  0-62  FP04  position  error  in  meters  for  case  4 
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Figure  0-63  Total  state  vector  error  in  meters  for  case  4 


The  50-run  average  position  errors  for  the  MN  location,  FP  04  location  and  total  state 
are  shown  in  Figure  0-61,  Figure  0-62,  Figure  0-63,  for  LML-EKFSLAM,  IMM-EKFSLAM,  LML- 
FastSLAM  and  IMM-FastSLAM.  It  must  be  noted  PEB  is  obtained  from  BFIM  which  is  calculated 
for  a  deterministic  mode  trajectory  as  already  outlined.  Averaged  MN  position  error  obtained  by 
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all  proposed  methods  varies  very  close  to  this  optimistic  PEB.  Unlike  expectation  from  the  single 
MC  simulation  results,  EKF-based  methods  have  better  performance  than  PF-based  methods. 
Specially,  in  terms  of  FP  location  error  and  overall  state  error,  EKF-based  methods  are  far 
superior  to  IMM-FastSLAM  and  LML  FastSLAM. 

Figure  0-64, Figure  0-65  and  Figure  0-66  illustrate  averaged  NEES  test  results  obtained 
for  all  four  proposed  methods  to  evaluate  the  consistency  of  solutions.  MN  location  estimates 
are  consistent  more  than  60%  of  trajectory  in  PF-based  methods.  However,  EKF  -based 
methods  have  overestimated  the  uncertainty  while  they  have  shown  better  performance  in 
terms  of  position  error.  In  contrast,  PF-based  methods  are,  as  explained  earlier,  underestimates 
the  FP  locations  uncertainty  that  has  a  dominant  effect  on  overall  estate  uncertainty. 
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Figure  0-64  MN  location  NEES  for  case  4 
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Figure  0-65  FP04  location  NEES  for  case  4 
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Figure  0-66  Total  state  vector  NEES  for  case  4 
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Discussion  of  results 


In  this  chapter,  novel  Bayesian  approaches  based  on  PF  and  EKF  has  been 
successfully  implemented  and  analysed  for  OWLS  according  to  four  system  models  defined  in 
Chapter  3.  Range  measurements  were  used  by  the  Bayesian  localization  framework  to  jointly 
estimate  the  MN  location,  stationary  FP  locations,  range  bias  and  ANs  sight  conditions.  The 
simulation  results  of  each  proposed  methods  were  evaluated  based  on  localization  error  on  a 
single  MC  run  ,  PEB  obtained  from  the  BFIM,  and  averaged  NEES  test  for  filter  consistency  and 
accuracy. 

In  the  first  case  model,  EKFSLAM  and  FastSLAM  were  applied  to  jointly  estimate  the 
MN  locations  and  FP  locations.  EKFSLAM  showed  comparable  localization  performances  to 
FastSLAM  method  in  terms  of  MN/FP  position  error;  however  FastSLAM  performance  were  not 
consistent  since  it  underestimates  the  FP  locations  uncertainty. 

In  the  second  model,  FastSLAM  and  EKFSLAM  were  extended  to  jointly  track  the  MN 
and  FPs  location  variables  as  well  as  single  range  offset.  For  the  same  number  of  particles  as 
previous  case,  FastSLAM  performance  degraded  in  terms  of  consistency  due  to  loss  of  particle 
diversity.  However,  both  FastSLAM  and  EKFSLAM  performances  outperformed  in  comparison 
with  the  first  case  solution  when  the  range  offset  is  ignored. 

In  the  third  model,  EKFSLAM  and  FastSLAM  were  generalized  to  jointly  track  location 
variables  and  range  offsets  from  ANs.  To  limit  the  dimension  of  state  space  which  particles 
represent,  the  range  offsets  were  tracked  jointly  with  FP  location  in  the  filtering  stage.  Unlike 
previous  methods,  EKFSLAM  consistency  degraded  where  it  overestimated  the  FP  locations 
and  overall  state  uncertainty,  however,  its  performance  was  still  superior  to  modified  FastSLAM. 

Finally,  the  last  system  model  were  considered  where  the  MN  receives  measurements 
from  unknown  sources  in  a  non-stationary  LOS/NLOS  environment.  The  system  model  was 
defined  by  JMNLS  where  sight  conditions  were  modeled  by  a  binary  state.  EKFSLAM  and 
FastSLAM  were  modified  based  on  two  sight  state  process  assumptions.  In  first  assumption,  it 
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was  considered  that  sight  state  process  evolves  based  on  Markov  chain;  in  this  case,  two  novel 
Bayesian  approaches  using  IMM  estimator  were  proposed:  IMM-EKFSLAM  and  IMM- 
FastSLAM.  In  second  assumption,  sight  states  at  each  time  step  are  independent  where 
EKFSLAM  and  FastSLAM  were  modified  to  detect  the  sight  state  using  LML  criteria,  called 
LML-EKFSLAM,  LML-FastSLAM.  All  four  algorithms  showed  better  performance  in  non¬ 
stationary  mixed  NLOS/LOS  propagation  in  an  opportunistic  unsynchronized  reception  in 
comparison  with  original  FastSLAM  and  EKFSLAM.  In  terms  of  consistency  and  position  error, 
proposed  EKF-based  methods  outperforms  PF-based  methods.  However,  the  LML-EKFSLAM 
provided  better  localization  error  in  terms  of  FP  locations  and  MN  position  error;  while  IMM- 
EKFSLAM  and  IMM-FastSLAM  had  better  performance  in  terms  of  sight  state  tracking. 
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